{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the basic libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython import display\n",
    "import os\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the environment for running the experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paths for loading/storing results.\n",
    "name = 'results/MixedGaussian_MINEE_MINE' # filename\n",
    "chkpt_name = name+'.pt'              # checkpoint\n",
    "fig_name = name+'.pdf'               # output figure\n",
    "\n",
    "# use GPU if available\n",
    "if torch.cuda.is_available(): \n",
    "    torch.set_default_tensor_type(torch.cuda.FloatTensor)\n",
    "else:\n",
    "    torch.set_default_tensor_type(torch.FloatTensor)\n",
    "\n",
    "# initialize random seed\n",
    "np.random.seed(0)\n",
    "torch.manual_seed(0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data using the mixed gaussian model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from data.mix_gaussian import MixedGaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_size = 400   # sample size\n",
    "rho = 0.9           # model parameter\n",
    "\n",
    "rep = 1             # number of repeated runs\n",
    "d = 1               # number of dimensions for X (and Y)\n",
    "\n",
    "X = np.zeros((rep,sample_size,d))\n",
    "Y = np.zeros((rep,sample_size,d))\n",
    "mg = MixedGaussian(sample_size=sample_size,rho1=rho,rho2=-rho)\n",
    "for i in range(rep):\n",
    "    for j in range(d):\n",
    "        data = mg.data\n",
    "        X[i,:,j] = data[:,0]\n",
    "        Y[i,:,j] = data[:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A plot of the first dimension of $Y$ against that of $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X/wZXV93/HnG5AN8iNKoCEoXzBoTCzDqN36Y2qijSRCJElDY2riiK0Td4hJJ2bsahXiT8wodBgz0TpuK9E1/mzF6AgGsG1CbRbr4hBcRFvWuAvCgrqSBaTg1333j3s+y/mePb9/3M85574eMwx7v/fccz7n3Hs+7/P5be6OiIjIEbETICIi46CAICIigAKCiIgkFBBERARQQBARkYQCgoiIAAoIEoGZ/bWZ/e6SjvV7ZnaPmT1gZj9RY/tvmdk5y0jbGJjZB83s0tjpkHFQQJBBJBnrQ0lGfI+Z/bmZHddwH2eYmZvZUS3T8BjgCuCX3f04d/9em/2U7N/N7Ml97lMkJgUEGdKvuvtxwDOBfwpcsuTj/yTwY8CtSz6uyCQpIMjg3P3bwOeBs7LvmdkRZnaJme0xs3vNbLuZ/Xjy9g3J/+9LShrPzfn8JjN7t5ndlfz37uRvPwN8I/X5/56XNjN7eXLs75nZxZn3nmVmO8zsPjO728zeY2ZHJ++FtP1dkrZ/ZWaPN7PPmdl3zOz7yb+fWHRdzOz1ZvZtM7vfzL5hZi+sOm7yvpvZq83s/yaffbuZnZl85oCZfTKVzheY2Z1m9kYz+25ScntZSZrON7Obk2P/rZmdXZVemRF313/6r/f/gG8B5yT/Po3FU/rbk9d/Dfxu8u9XArcDPw0cB1wFfDh57wzAgaNKjvM24EbgHwEnA3+bOk7p54GnAQ8AvwBsYlG9tJ5K9z8BngMclezrNuA1qc878OTU658A/iXwWOB44L8Af1lw7KcCdwCnptJ6ZoPjfhY4AfjHwMPAf0uu4Y8DXwNekWz7guScrkjO8fnAg8BTk/c/CFya/PuZwL3As4EjgVck3+OmsvTqv/n8pxKCDOkvzew+4IvA3wB/krPNy4Ar3P2b7v4A8AbgpQ3aDV4GvM3d73X37wBvBV5e87O/CXzO3W9w94eBPwYOhjfd/SZ3v9Hd1939W8D7WWSoudz9e+7+KXf/gbvfD7yjZPsfschon2Zmj3H3b7n77gbHfZe7H3D3W4FdwHXJNfwHFqWxZ2S2/2N3f9jd/wa4GvitnDS9Cni/u3/J3X/k7h9iEWyeU5ZemQ8FBBnSv3D3x7n76e7+and/KGebU4E9qdd7WDwZ/2TNY+R9/tQGn70jvHD3B4FDDc9m9jNJtc8+MzvAIqCdVLQzM3usmb0/qYI6wKLK63FmdmR2W3e/HXgN8BbgXjP7uJmd2uC496T+/VDO63QD/veTcwuKrtHpwGuT6qL7kmB+GotSQWF6ZT4UECS2u1hkRMEaiyqOe1hUjbT5/F01j303iwwPWGToLKp9gvcBXwee4u4nAG8ErGR/r2VRtfLsZPtfCLvO29jdP+ruz0vS78C7Wh63yuPN7NjU66JrdAfwjiSIh/8e6+4fq0ivzIQCgsT2MeCPzOxJSbfUPwE+4e7rwHdYVOH8dMXnLzGzk83sJOBNwF/UPPZ/Bc43s+cljbBvY+M9cTxwAHjAzH4W+L3M5+/JpO14Fk/n95nZicCbiw5sZk81s180s03A/0s+96Oax23jrWZ2tJn9PHA+i/aNrP8EXGRmz7aFY83sxWZ2fEV6ZSYUECS2K4EPs6he+XsWmc2/BXD3H7Coh/9fSRXGc3I+fymwE7gF+CrwleRvlZL6998HPsqitPB94M7UJv8O+B3gfhaZ5Scyu3gL8KEkbb8FvBs4Bvgui4buvyo5/Cbgncm2+1g0ir+x5nGb2sfi3O4CPgJc5O5fz27k7jtZtCO8J9n+duBf10ivzIS5a4EckbkysxcAf+Huhd1fRQKVEEREBIgYEMzsx8zsf5vZ35nZrWb21lhpERGRiFVGZmbAse7+gC3mnPki8IfufmOUBImIrLhWk4b1wReR6IHk5WOS/9SgISISSbSAAJAM2LkJeDLwXnf/Utn2J510kp9xxhnLSJqIyGzcdNNN33X3k6u2ixoQ3P1HwNPN7HHAp83sLHffld7GzLYAWwDW1tbYuXNnhJSKiEyXme2p3mokvYzc/T4WE56dm/PeNnff7O6bTz65MsCJiEhLMXsZnZyUDDCzY4BzWAzXFxGRCGJWGf0Ui1GeR7IITJ90989FTI+IyEqL2cvoFg6foldERCIZRRuCiIjEp4AgIiKAAoKIyOhs3b6Drdt3LP24CggiIgJEHpgmIiKPCqWCW/bs3/D68gufu5Tjq4QgIhJJrKqhIioh1LTsSC0iqyfkL7HyGwUEEZEli101VEQBocJYvzgRma9Y+YsCgojIksWuGiqigFBhrF+ciEjfFBBERCIZ2wOmAkJNY/viRORRKsH3Q+MQREQEUAlBRCZMvQD7pRKCiIgAKiGIyISpF2C/VEIQERFAJQQRmQGVDPqhEoKIiAAKCCIiklBAEBERQAFBREQSCggiIgIoIIiISEIBQUSiG9vawqtKAUFERAANTBORiDQ53biohCAiIkDEEoKZnQZsB04BDgLb3P1PY6VHRJZPk9ONS8wqo3Xgte7+FTM7HrjJzK53969FTJOISFQxg2O0gODudwN3J/++38xuA54ADBIQ5voEMtfzktWi3+84jKJR2czOAJ4BfCnnvS3AFoC1tbWlpks2UvCRuYv5Gx9DA3v0gGBmxwGfAl7j7gey77v7NmAbwObNm73p/sdwkYcw1/MSkXiiBgQzewyLYPARd78qZlqk2JDBR4FMxmAMD1hjaGCP2cvIgA8At7n7FUMdZwwXeQjLOq+t23ewe98BzjzlhEH2LyLjEbOE8M+AlwNfNbObk7+90d2viZgmyXHmKSdw+YXPHaRkoCovGUqT39SYHhxjHjtmL6MvAras4801oxmyZAAbM2yVFETmLXqjsgxjiCedUFLow5ieyGReupQ+V/13qIAguZRhy1Tt3ndYZ0WpSQFhILEy0qnVzY81XTINeb/vbLWmfmP1KSCsmKYBQjeTjEnZ7zf7MHTspu7Z29gfqPqmgNCz2E/oquqRZRpTSTjb6UEdIJpTQFgRsQOVSBd1nv777B69qveLAkLPxvKEHuuJbe43jCzEzjDVTjAMBYQVsYxAtXvfAbZu36Gbc+LGGNyLfr956zD3ke6xPNgtmwLCQFblB5R9UlRQWA3LzDDLjqHfWb8UEFbMUCWD4MGH1xUUJip2NVAd2bRUpbHrOYzp3JdBAUE6CY14u/cd4MGH1wH17lglyygZjDlAzY0CgnSWDgoKBtPVZzXQMmbhheJgoWDSjgKC9CLd3U+kD6vasBuTAoL0Tk9l0zaF/vtVwULBpB0FBBEZTB8TzSkzXx5zb7xMcTSbN2/2nTt3xk6GVIjxVKYnwXG54LJrgUc7GOh7icvMbnL3zVXbqYQgIp1l249Cj7OqqiMF8nFRQJDexSgZqN2if02vZageUk+z6VJAkF7VmZ5YmfVwhrrGVU/42Unnzj79xA3bVX2uS7r1u+qPAoJMmnqT9K/petp5I9VBJYUpUqOy9CKbiYQnxPT4hPRTZNX6zE0z+LrbzzVwlF3/rvs6dtNRh9oE8vabHane5Nh9lAz6OOe5U6OyrBRlAv1Jl7pCySBkukXbp7fVdzFdKiFIr8rqmus8RQ711FdU133V617Uab9j02cJKL2von8Pdeyq9CzzuHOgEoKI1FI12rcuZcjTpxKCLFXsuv4wYCpbUglWMVNr+p0su85ebQXdqYQgvVPRfF6WNYZDv5vpUECQpaqbKQyVeYQ2g+zI2qkPbFtGuofo4ltnX+pavDwKCFJpCqOBx5imIfVxvkNntG1+N6v2PY5N1IBgZlcC5wP3uvtZMdMi/YndTlDHXHqqZNtElllS6CIbLC647NrKgWxNjzvV7zSm2CWEDwLvAbZHToeUWHaRvclxupZepraoT/Z8j7Du+xzq+2zyu3nokfUN3ZKVmccRNSC4+w1mdkbMNEh/6mbOY6yCmmrGczDpJBjGVQxxHkN8P2E8QxgFfbDHzo5j/H1NRewSQiUz2wJsAVhbW4ucmtXQV7/0tsdtciPXeQrNey87UC5WptH0uGG7bFXRmNU9tzNPOYFde/dzzNFHKfOOZPQBwd23AdtgMQ4hcnKkRN0qgqLtyqpvqvbZx8pcU5Ktby+6Ln3OFVRntHITeb+DEOi6UK+k9kYfEGR5ll3UzmbiXW7kvAbJovMJwtP1kNUtRbpe6zpBtGlaYmacu/cdYOv2HdFLbKtOAUF616SKoEze5GpFpYqp1Rf3VaKpKhl0uS55AXrr9h1s3b6j15XQ0iWOPo39NzBGsbudfgx4AXCSmd0JvNndPxAzTatsWUXtqsyq7Om3SUaabrjMm4VzmcEje6y6VT7LSNMQwbRpwFM1zzjE7mX02zGPL8vRtitoesrlkKlnhUw/3f2yz8yk7wyqaNbVvo/XZwab/mzVfsPfVfUzTaoyksMMffO2yayK+qnD4mk0dFsMddFBOqhs3b4jN3MbUlUAGHpVsbIn9SGeyrPHa1tSkDgUEGQwbaskiuqu0/uCRzObdNfLvpZwHKo6paiqaFnH60tRusLxwnloGc1pUUCQaPpo4ISNQQIWmVAICkNOldxlnESbRtShRnAPcW366rmlKqflUkCQwXStkmhSkkiXCvrotTJ0I2fRoL+5ZIAqGUyTAoJMSlVVRWhDyC6m0lVZyaRJSaHKrr2Hr1081AjuvGO0DUTZNHa97lPtTjx1CggyuGU1UkP/A7X6HABWxzFH93NLdklvm8y3z5Hi6RKfLJcCgsxKX1UvRb2V+p6+IcjOTRReX/W6F5WeU91pQorkPYmHwYBN9D2uIl3iyxtDIsNQQJBZ6ePJOGSOffRWWrYmVS15mX/o2nvLnv2H1igoy4yrutV2SXs6TdkuwzIMBQSZpS4lg+DBh9cPZXB5vZX6zKDC0p7pkgHkV1+l/95kgZk86Ub47DQhDz1SPZNqtmpniOA5pYA8dQoIMgt9zd2TnRYbaFWF0jYNXdWpMssLJg89ss5ZayduyOAPevnTedhP+tr0OftpH/uUZhQQZNLqVBE1nWjtgsuu5QiDs9ZO3PD3ZciWDMrmewrzND348PqGEdxVac1rsA1rEJx36dUb/h5GiGelA+cte/Z3qiqS8dC3KIdMuWtfn/34Q7XJrr37D02JUbXfrt1S29i970BhtU7ZsdINtvBoQ/Z5l1592Mpl6fPPjvvILs5TpxTVZxdd6Z8CgoxeWe+aorUOyrapO2V03bR0UVQ9UifAlfXsqXPNuqzHnG5rKJpNVqZHAUFmNQioa5qzdehHGBuWdExn4HkNsemeOUNdx7569qSrxODR6rKHHlnfsFZzNrPPlhTq9ETKW9ei7LNT/g1OmQKCRNP0ib1swFjZHD1NJtUrqg7J/r3O4Kmy94vWRM5rK2iiyTXLk54HKryuqn7KC37K0KdJAUFW9ibOO990UAjy6sbTwWH3vgMcYYsSRd703H0LaQyll6q6+6p+/EVP/9n3yj5XdtzsuhZBaJCO0f4i+RQQZOnq3vRNSgFl2mYmIXPftXc/xxx91KHXVXXv2cbavAwvu55zyNTrpDVdrRMy1aLjlE30V3b9616zvC6s6fPTFBTTooAgh6xKXW7Tp9CDvnGQVvapPBu4sr1/+p7nJ13HnxWeykMa0k/oy5qPKS1UQdVdwnRVS6tjoYAgS9embr/J9l3ltSWEQVvpjK0ogy3LsLucS0hTet/paqN0aSDbFtBXGuruq04JRcZHAUFylT1Fj+XprW066lZFZXvCVI2gDdtkg8FDj6z3Notpkbz6+iB2t9A6bRdN3pPhKCBING3rqZs+cXYJHOljVX2+qOdRXsNvmwwvPYo6BJ1019G6VVNNrkfTUlzRa5kGBQTJVVT0z6uTXvbNX9Zo28QQDdLpvvawqNIJ01Gk1ZlrqGgwWnZ1uPTnsvMuKWOWJhQQZPSK6uzLRhSnt286UrlNJpoNnGHls+zgr7qy6zEE6cFoITBkz7No2u4m10PdP1eTAoKUGmMPkLGko0y6oTetTttMnfUYysYeNOnCKpJm7gXdIUZo8+bNvnPnztjJWHljyYjrpqPpiOi8tQ+q9p2dRiJUZxWtLZx3rFCyyOuhlN1PXtrK1lLIS3OfbQh1jeW3s2rM7CZ331y1nUoI0thYbubY6cirsy8SnvbLMuu8uYSK9lU28rgOrUImeQpLCGZ2DfBqd//WUlNUQiUEGVLTp9eykcHZQJEtEQTZ6q/0TKTprqrphumiPv9VpZRsA3STc+2qSylMuuujhPBB4Doz+xBwmbv/sK/EiXQ1ZIbWtIopveJYVrYraGhs/vwlLy7dLm/cQtfut3X20ef1VPXQ9BQGBHf/pJldDbwJ2GlmHwYOpt6/ouvBzexc4E+BI4H/7O7v7LpPkSp15k5qU6VSNDgs3RW1aIBa6E4aBpEFedNFZ9UdLRxzxPAUOgJIdRvCD4EHgU3A8aQCQldmdiTwXuCXgDuBL5vZZ939a30dQ+ZnyO6Qeb18skEhfbzsRG7Z5SbTo5zD0pRFs6Gmg0hRxt/HOgtVXXizo7PbUJfV6SoMCMnT+xXAZ4FnuvsPej72s4Db3f2byfE+Dvw6oIAgg6iTUaUz9AcfXq8sKaQnsUsvKgMbB46FkkF2/YPsMcM+d+87wFWve9FhGXedkkL6ddFgwlgUFMatrIRwMfASd791oGM/Abgj9fpO4NkDHUtmYhlVD2Hheni0T39eMAnVPNleQeF1eq7/ot5FYfbS9HQUYR9ZRSOU28jrirpr7/5DE/pl1ylos2+VDKanrA3h5wc+dt6s8od1eTKzLcAWgLW1tYGTJHNWJ6NKL0BfZ2bTPOm1E/JkSx3Z0kN6Yfvs0pNAo8xambM0EXMcwp3AaanXTwTuym7k7tuAbbDodrqcpMnYDdEbpqhraPp42cCQzfirJrELbQ+79u7fUJVz9uknHtYYHUoEeZPWpaet6HotwjmFBXdClVffpRAZv5gB4cvAU8zsScC3gZcCvxMxPbIiYmVUIePNW00s/e9jNx3FQ4+sbyitpKuxsj2R6lDmLHVECwjuvm5mfwBcy6Lb6ZUDtleIFGqzqEv4TLanUfb9tLwn/fQCN6GEkG6DSDdQZz/fVy8eVStJEHXqCne/BrgmZhpk9RT114fDxxJUZZJ11zkI+0ln/NnlJfPaDKqqpET6pLmMRBJlJYOqLqt1Gp1Dxh9KFUWN2kVzGmUbufteCU0lA1FAkMlrOutpyNTPu/TqDT2CiuYmSr9XpKxkkF6bOV1SKNtH2UI4Y6OqpvlQQBCpoc969qoMPhtE0t1Us0EtbD9UZqzMfrUoIMhkNZ0iIWSoocfOQT+8Hj+9bdU+606Cl+4dVLZ9mbrrJS+TpqiYHwUEkQaaTnbXVFEQya6CNnTmu8xZUWU8FBBkstpkjnW7mDZZ6jK7UllQ1AOpSXrTXVCLjlNX30EkOztsn/uWOBQQZNTqZGJjq04pWrCmjmwmm522IrtdnbTUXdWtKB1pqh6aN62pLKNWJ+NpmznVWcWrbN/ZQWnZKS+KpsJIB4qiBuayabaz6zAXpS/sZ9fe/Rt6U7VZM7ruOWkltHHSmsoyaXUaLMfaqFm1YE2Zvs4h21OpbLK9Mnmzoqb/Hns6bemXAoI0MpZMtw916vTLzrNoSus8YZvd+w4cmuo6PTah7LN51T1NR1TDo9Nk9PndadqLeVFAkFFqkll3qVZqW8ee3mddYe2Esn1VZfB120vS1yZUGeWdY9OMXBn+vCkgSC1jrJ7pq7oib92DJudV1eCdfpIPsmMSsucSgkf6s2EW1PR0Ftl2iLalna6qqpZkGhQQZNSadCXNU9WfPruWMFRPIJfX7TQvUw7v79p7eC+hhx5ZP7QwTvocshPtZT+TDgZVpYXsOgdlXWebTv+hjH6eFBCkliHriuuO+C3KcMvmIKorvRRml+Uj8xxz9MYn+7LeRSGTTzcCpwenhXM+IllvMNuOkQ1wy54dtUmgUXAZHwUEmbWqQBaqZLK9cMoy0+wUGNnG4aLBb3UmqwvHzBt7AItunXVLMulptbPnXtUTqqgKaExVhtI/BQRpZIiSQdWU0lXvh0yvaa+btPQi9+kn+C5jHEJjbjqNbfaVzvjTJYi8ksyyMuqi6qomc0AVVblJPAoIshLKBpbBIhiEoFAnc6qT8WV79tSpn08HETi8CqnuqOd0usqOWzf4atzBatBIZYmuaRtCmfR8P1Wfy44ADvXyn7/kxbVGMeelLTsgDMrbDLL7KqqqyhshnDdLa1G6qhSdbxD2UzZquq509VmX/Uh9GqksK2fr9h2H9d5JvwcbM5zQINtl0rg+MrBsELllz/5Dg8iKBr8VPam3qesva5TPHqeqjUOmTQFBoqtbPVMmZKqhi+UFl127obG3aY+btnXyeZPTtQkaZQPJ0nMTde0RVXRtmrTfNDlWep+hgVztB+OhgCCTl1dNkx4RHN7LZp7h322nk26TznDstGwQgcPbMerW4Zdl2GXVW9lrk25fkdWhgCCzkO6jH+rs4dF69rZVHG2fXNMBpw95k9V1mZsoL4imSwphJbmi9o82x9S8R+OngCCjVbfapaiaZuv2HbndR7t2UW1zHlCvXj8Er65VQXklg7xBa3lBNH1tjrDxrTchwzkidgJEipRNBlckHTyqulwOJR1wbtmzn937DnTOVC+/8Llc9boXcfbpJ3LspqM4+/QT+fwlL+bMU05oVRIJ1ynsL33d0mkNYzOqjnHBZdfWrmbKfi/hekl8KiHI6OT1uukyeKmqT/3QQaNoKc2sodJVtd9sEM1e/yaN8dmMfehSl6qd+qWAIJ31fXPmlQzaLvASQ9W0EH3tu49qr6JG57rVddlxCWEiv7PWTiz8TPZ4mg5jPBQQpJEu3Q3rfibb172q8bQoYwmqRt8uS9Ouq7GOH7ZtE8BCFdNQmbyCyLAUEKS1oW7O8PnzLr2ag75x+ubQB39ZXUW7qBrVnKfuNeyzeqnLYLbwPZx36dUAG6bnrlJVklImv3wKCFJLWU+VJp+B+jd6GHFcVV1UlbEMMS3C1IJTV1XfXfiu6raXtKWgMawoAcHMXgK8Bfg54FnurgmKJqjrjKBV0tM27Nq7/7CRyOltitSdKjoY4lzqBMa2I4P7SGeTTLYondlpQLqMCld1UDyxSgi7gAuA90c6/qwNcSO1eTJb9tNcUcZSZ4bQuvXlITiFqpEHH17nvEuvnkRJoeh7KDv3ptcyBIJlfdfSrygBwd1vAzCzGIeXni3z5m9aMoBFph3WNi5Ka3b2zbIg1qSePL2Psn1WVXst46m5zj6LqoT6SKeqg+IbfRuCmW0BtgCsra1FTs24LSPzmMJN2qQvfbhGdbu1pvcdZladSsmgaS+s9L/LutDWaU+SaRgsIJjZF4BTct662N0/U3c/7r4N2AaL9RB6Sp5MVN3Mt0lf+rRQJZK3fTZjBRqNpG5SzZZ9Heupuag6qWhwWx/pnMJDx1wNFhDc/Zyh9i35YmceWX2VWLoMumo6O2gT2UVkxqrqd1Hn+tYJkCopTN/oq4ykm7EEh1i6jNwtei/GNY1VMmgazLW2wbTF6nb6G8CfAScDV5vZze4+7orYCRnLDVknA63z3lgaVKeu6fWvu79Vf+iYk1i9jD4NfDrGsVeF+nQPZxWuoTL71aQqIzlkqJu/7dO/MqXh9PnAMMbvRb+ZdhQQZkqZ6bhN5XvpO31TOe9VpYAgUaqXmgQsZR79m+sDg6pKu1FAmDndCOOyqhlWlxlV535txkQBQSbZnVKZRT/mdv3mWvJZFgUEmayq+YnGaFUzrHCeYS6qNjOqrsq1ikkBQQ6Zwg2XzSymGBTmYswZ9RjTNAUKCDI5TWcyHaMppbUP2UkE28z82vaYq3atu1BAkElpMpOpDENVOvOlgCBL1ddsmG1mMpV40k/9oYRX9b11LRkoYDWngCCtxbzR6sxkKsMYqmFcGXd8CggjMfebYYintrleq7lKtyPcsmf/4FOlzP2eGoICggDNbh4VyQX6Lxno9xSfAkJkq3Iz6KlN6izHOcTxpD4FhAnqM1NtE5CUuUuf+v496XfZngJCZMvIXMd0g4whDRLXGH8DY7pHYlJAmJAhG2bb7GvVbx5pp+i31lfJYO7Vr0NSQBiJIUsGukFE8uke2UgBYUKGrF5a1RtAlmfozLeP+yM9LcoqUkCYMTX+yhCW/Xsa8nhD93SaGgWECVLGLlO0rAeULiWDOpPvzZkCwoDG8qOKfXyZh2XXty/zeGGCxHCsVaWAMFFjCTYiTfU5fmZM7Q9zoIAwAPVckDlqm2m2/f2PPZMea7q6UECYGAUbWWXL6qm0qhQQBjD2J5sp0LUbr6Ylg66Z99h+A3N+KFNAmBgFG5krzaMVnwLCgPRjbW7OT1+rZq6Z91zPCyIFBDO7HPhV4BFgN/Bv3P2+GGmZqjn9CGW1NXkImGMmPCaxSgjXA29w93UzexfwBuD1kdIiIzLnp69VNdfvcI6/1SgBwd2vS728EfjNGOmQw83pxy3TUCdjHaIqUb/1w42hDeGVwCeK3jSzLcAWgLW1tWWlSSLTTSpjN8f2LnP3YXZs9gXglJy3Lnb3zyTbXAxsBi7wGgnZvHmz79y5s9+EdjSHHwEc/uM++/QTgemfl8xL1/tt6/Yd7N53gDNPOSH3t95lbfE690ys/MLMbnL3zVXbDVZCcPdzyt43s1cA5wMvrBMMJI65BDyRvi2rDWGZ92CsXkbnsmhEfr67/yBGGrqaW3Fx2Qugi7TRpWQAGyevO3bTUZx5ygmHSgZbt+8Y7H6eSn4Rqw3hPcAm4HozA7jR3S+KlBbJMZUfsEhsbafbDjOsFolxD8bqZfTkGMft0xDFxTFkuioZyBxV3a9DV/+k9x+CwRgfrsbQy0hGaI59rEViCsHgwYfXuWXP/sp7K8Y9qIDQUZ8lA1XPiAyr6p4a+p5L924aIwUEKaWgJNIn3zAyAAAEgklEQVSPtk/8y7wHFRBGQNUzIqtl974DbN2+Y3T3+hGxEyAiskouv/C5lT2MYhlspPIQxjhSWUSkrlgzAtQdqawSgoiIAGpDEBFZmrG3F6qEICIigEoIkzTWpwsRqWes965KCCIiAqiEMCl1RjSr9CAibamEICIigMYhTFJZyUArnolIlsYhiIhIIyohzIzaEEQkSyUEERFpRL2MZkYlAxFpSyUEEREBFBBERCShgCAiIoACgoiIJBQQREQEUEAQEZGEAoKIiAATG6lsZt8B9jT4yEnAdwdKzhjo/KZrzucGOr+xOd3dT67aaFIBoSkz21lnuPZU6fyma87nBjq/qVKVkYiIAAoIIiKSmHtA2BY7AQPT+U3XnM8NdH6TNOs2BBERqW/uJQQREalJAUFERIAVCAhm9nYzu8XMbjaz68zs1Nhp6pOZXW5mX0/O8dNm9rjYaeqLmb3EzG41s4NmNpsufmZ2rpl9w8xuN7N/Hzs9fTKzK83sXjPbFTstQzCz08zsf5jZbclv8w9jp6lPsw8IwOXufra7Px34HPCm2Anq2fXAWe5+NvB/gDdETk+fdgEXADfETkhfzOxI4L3AecDTgN82s6fFTVWvPgicGzsRA1oHXuvuPwc8B/j9OX1/sw8I7n4g9fJYYFat6O5+nbuvJy9vBJ4YMz19cvfb3P0bsdPRs2cBt7v7N939EeDjwK9HTlNv3P0GYH/sdAzF3e92968k/74fuA14QtxU9WclltA0s3cAFwL/APzzyMkZ0iuBT8ROhJR6AnBH6vWdwLMjpUU6MLMzgGcAX4qbkv7MIiCY2ReAU3LeutjdP+PuFwMXm9kbgD8A3rzUBHZUdX7JNhezKM5+ZJlp66rOuc2M5fxtVqXWVWBmxwGfAl6TqYWYtFkEBHc/p+amHwWuZmIBoer8zOwVwPnAC31iA0safHdzcSdwWur1E4G7IqVFWjCzx7AIBh9x96tip6dPs29DMLOnpF7+GvD1WGkZgpmdC7we+DV3/0Hs9EilLwNPMbMnmdnRwEuBz0ZOk9RkZgZ8ALjN3a+InZ6+zX6kspl9CngqcJDF1NkXufu346aqP2Z2O7AJ+F7ypxvd/aKISeqNmf0G8GfAycB9wM3u/qK4qerOzH4FeDdwJHClu78jcpJ6Y2YfA17AYnroe4A3u/sHoiaqR2b2POB/Al9lkacAvNHdr4mXqv7MPiCIiEg9s68yEhGRehQQREQEUEAQEZGEAoKIiAAKCCIiklBAEGkpmfny783sxOT145PXp8dOm0gbCggiLbn7HcD7gHcmf3onsM3d98RLlUh7Gocg0kEyjcFNwJXAq4BnJLOYikzOLOYyEonF3X9oZluBvwJ+WcFApkxVRiLdnQfcDZwVOyEiXSggiHRgZk8HfonF6ll/ZGY/FTlJIq0pIIi0lMx8+T4Wc+LvBS4H/kPcVIm0p4Ag0t6rgL3ufn3y+j8CP2tmz4+YJpHW1MtIREQAlRBERCShgCAiIoACgoiIJBQQREQEUEAQEZGEAoKIiAAKCCIikvj/3yrkc5LSNHAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[0,:,0],Y[0,:,0],label=\"data\",marker=\"+\",color=\"steelblue\")\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.title('Plot of data samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the MINEE model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from model.minee_mine import MINEE_MINE "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_size = 100       # batch size of data sample\n",
    "ref_batch_factor = 10  # batch size expansion factor for reference sample\n",
    "lr = 1e-4              # learning rate\n",
    "\n",
    "minee_mine_list = []\n",
    "for i in range(rep):\n",
    "    minee_mine_list.append(MINEE_MINE(torch.Tensor(X[i]),torch.Tensor(Y[i]),\n",
    "                            batch_size=batch_size,ref_batch_factor=ref_batch_factor,lr=lr))\n",
    "dXY_list = np.zeros((rep,0))\n",
    "dX_list = np.zeros((rep,0))\n",
    "dY_list = np.zeros((rep,0))\n",
    "mine_dXY_list = np.zeros((rep,0))\n",
    "mi_list = np.zeros((rep,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Previous results loaded.\n"
     ]
    }
   ],
   "source": [
    "load_available = True # set to False to prevent loading previous results\n",
    "if load_available and os.path.exists(chkpt_name):\n",
    "    checkpoint = torch.load(\n",
    "        chkpt_name, map_location='cuda' if torch.cuda.is_available() else 'cpu')\n",
    "    dXY_list = checkpoint['dXY_list']\n",
    "    dX_list = checkpoint['dX_list']\n",
    "    dY_list = checkpoint['dY_list']\n",
    "    mine_dXY_list = checkpoint['mine_dXY_list']\n",
    "    mi_list = checkpoint['mi_list']\n",
    "    minee_mine_state_list = checkpoint['minee_mine_state_list']\n",
    "    for i in range(rep):\n",
    "        minee_mine_list[i].load_state_dict(minee_mine_state_list[i])\n",
    "    print('Previous results loaded.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuously train the model with MINEE. The following can be executed repeatedly and after loading previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4VFX6wPHvOyU9QIAAAUSUoiAKAirYseyqq2tZ29pdlSZFxQYESICAXUC6va76k7V3saKigKCCqBTpSQgQQkmbcn5/3GEyk0ySSUgyJHk/z8PD3HvPvfedycw7Z8499xwxxqCUUqphsUU6AKWUUjVPk7tSSjVAmtyVUqoB0uSulFINkCZ3pZRqgDS5K6VUA6TJvYERkS9F5NYIx3CUiCwXkb0iMiKM8mki8pLvcQcR2Sci9tqPtOESkWtF5JNIx6EiR5N7PSQiG0SkwJcEs0XkWRFJqOIxOoqIERFHLYR4L/ClMSbRGDOjKjsaYzYZYxKMMZ5aiKtBCvW3NMa8bIz5Wy2dL+IVCFU5Te7110XGmASgN3ACkBrheAIdDqyKdBCBaulLTKlDlib3es4YsxX4EOhRepuI2EQkVUQ2ish2EXlBRJr6Nn/t+3+37xdAfxHpLCJfiUieiOwQkdfKO6+I/FNEVonIbl9Nrptv/efAAGCm77hdQ+x7hO88e0XkU6BlwDZ/LVRErhaRpaX2vVNE3vE9jhaRR0Rkk+8XzFwRifVtO1NEtojIfSKSBTzrW3+viGSKyDYRudV3rs5VON4o32uZKSI3B8QVKyKP+l7rPBFZFLBvPxH5zvda/SwiZ1bwurYVkQUikiMifwU2a4nIiSKyVET2+OJ7rIK/5U0isihgXyMiQ0Vkje91nyQinUTke9/xXheRKF/ZJBF5zxdDru9xe9+2DOC0gL/vTN/6o0XkUxHZJSJ/iMiVAee+QER+8513q4jcXd7zVzXIGKP/6tk/YANwju/xYVi15Em+5S+BW32P/wOsBY4EEoD/AS/6tnUEDOAIOO5/gbFYX/oxwKnlnL8rsB84F3BiNcOsBaJKx1DO/t8DjwHRwOnAXuCl0nEBcb5tXQL2XQJc7Xs8DXgHaA4kAu8CU33bzgTcwIO+88QC5wFZwDG+Y7/oO1fnKhxvou85XwDkA0m+7bN8z7sdYAdO9p23HbDTV97me812AskhXhcbsAwYD0T5/m7rgb8HvG7X+x4nAP0q+FveBCwKWDa+59bE9/yLgIW+czQFfgNu9JVtAfzL9xolAv8HvBVwrKC/LxAPbAZu9v3degM7gGN82zOB03yPk4Dekf4MNYZ/EQ9A/1Xjj2Yl933AbmAjMBuI9W3zf/B8H96hAfsdBbh8H8BQCeEFYD7QvpLzjwNeD1i2AVuBM0vHEGLfDr4kGR+w7hVCJHff8kvAeN/jLljJPg4QrC+YTgHH6Q/85Xt8JlAMxARsfwZfsvYtd/adq3OYxyso9XptB/r5nn8B0DPE870P3xdqwLqPDyTSUutPAjaVWjcaeNb3+GsgHWhZqkyov+VNlE3upwQsLwPuC1h+FJhWzt+sF5AbsBz09wWuAr4ptc88YILv8SZgENAk0p+dxvRPm2Xqr0uMMc2MMYcbY4YaYwpClGmLlfwP2IiV2FuXc8x7sZLcj74ml/+UUy7ouMYYL1bNrV0YcbfFShT7S8VVnleAf/seX4NVg8wHkrGS/DJfc8du4CPf+gNyjDGFpc69OWA58HE4x9tpjHEHLOdj1aBbYv3SWRci/sOBKw4c03fcU4GUcsq2LVV2DCV/r1uwfjX9LiJLROTCEMeoSHbA44IQywkAIhInIvN8TUx7sL5Umkn5PZgOB04qFfe1QBvf9n9h/XLZ6GuO61/FuFU16EWmhm0b1gfvgAO15mxCJGJjTBZwG4CInAp8JiJfG2PWhjjusQcWRESwmoe2hhFTJpAkIvEBCb4DVs0ylE+AliLSCyvJ3+lbvwMrIR1jrOsOoZQ+ZibQPmD5sIDH4RyvPDuAQqAT8HOpbZuxau63hXGczVi/FLqE2miMWQP8W0RswGXAGyLSgvJfu+oahfUr7yRjTJbvtV+O9cVPiPNtBr4yxpxbTtxLgItFxAkMA14n+LVXtUBr7g3bf4E7xbqAmQBMAV7z1T5zAC9WmysAInLFgQtnQC7WhzhUl8TXgX+IyNm+D+worDbc7yoLyBizEVgKpItIlO9L5KIKyruBN4CHsdrCP/Wt9wJPAo+LSCtf/O1E5O8VnP514GYR6SYicVht2wfOU53jBe77DPCY74Ko3XdRMxqrWekiEfm7b32M7+Js+xCH+hHYI9ZF4Fhf+R4icoIvnutEJNl3vt2+fTyE+FsepESsL7rdItIcmFBqe3apc70HdBWR60XE6ft3gu91jhKrz31TY4wL2EPo95SqYZrcG7ZnsC4afg38hVW7HA7ga9rIAL71/ZTuh9Wl8gcR2Yd18W2kMeav0gc1xvwBXAc8gVVrvQira2ZxmHFdg9W+vAsrcbxQSflXgHOA/yvVLHIf1oXcxb7mg8+wapwhGWM+BGYAX/j2+963qag6xyvlbuBXrAu+u7Au5NqMMZuBi7GaV3Kwarn3EOKzZ6y+/RdhtXH/hfXaPoV1wROsC8KrfH+f6VgXlgvL+VsejGlYF6B3AIuxmqcCTQcu9/WkmWGM2Qv8Dbga61ddFiUXsgGuBzb4XtPBWO8dVcvEGJ2sQzVOYnXfXAlEl/rSUKre05q7alRE5FJfU0ESVu3yXU3sqiHS5K4am0FYzSPrsNp+h0Q2HKVqhzbLKKVUA6Q1d6WUaoAi1s+9ZcuWpmPHjpE6vVJK1UvLli3bYYxJrqxcxJJ7x44dWbp0aeUFlVJK+YlIRXd0+2mzjFJKNUCa3JVSqgHS5K6UUg2QJnellGqANLkrpVQDVGlyF5FnxJpWbGU520VEZojIWhH5RUR613yYSimlqiKcmvtzWKPRled8rBlyugADgTkHH5ZSSqmDUWlyN8Z8jTWEaXkuBl4wlsVYM7aEmmVGKaUOee++NovnH70/0mEctJq4iakdwdOVbfGtyyxdUEQGYtXu6dChQw2cWimlatay1TlYsybWbzVxQVVCrAs5GpkxZr4xpq8xpm9ycqV3zyqlGqHcnTuZm3Y7uTt3RjqUsLmKCvjoycfwuA+d0aNrIrlvIXg+xPZYs7EopVSVvT59Clkk8/qM9EiHUqmNPy8iLS2NmZMfZfHWPbw/77EKyxcXF5Ofn18nsdVEs8w7wDAReRVr6rQ8Y0yZJhmllAqH2xsLNjDe6MoLV+LlWVNZk1NEW3IZmDa9zPZvPn6L5d8s57yrL6Nrj55lts9MHUdCMxc33f1AyOP/vMiaNjjPbtXY9++veKbJWbNmkZeXR1paWhWfSdWF0xXyv1hzTR4lIltE5BYRGSwig31FPgDWY809+SQwtNaiVUopnyVff8CjYyfxwmNjyi2zPycLgH2uFiG3r/j6J3Y5hC9eDz2N7w6HnQ37Ynh3weyQ26V0o3SoRuoAeXl5FReoQZXW3I0x/65kuwFur7GIlFL13sdvvcCapWtofWQSV/znrrD2mT9+JJm0oHllGdJn+cKP2Otsjm1nycXP+enDAYhJbEbB7gKMiQI77HF6eWXOw1wz5B6efXQMRbs9DJ70YIXHD6xd71i/ocz2J1LTKcARlEX3FQjPZjxGbmExl11yNh37nMAX739PbFws/Qb0Cut51ZSIDfmrlGq4/lr2Kzsc8djW7wh7n0xaYGxeXB4n4AVgTtrtiCeewZMeCusY24yvhr4HsCUEbVuTWQjAxj1Wwp+ZmkoBMUB4F0F3bd/MW3NncvWoMcTFN2WnwwCuoDJbvS7wusAO//fW5+x/933/Nokr8j8uKCggNjY2rPNWlyZ3pVS53nnjGTauXMKZl9zEsb1O8q/P3raNBfMmE9+qNXsz3XTt25m/XXJjwJ5W7duEUQufN2EEXjwYm9WDbo/TW3IeksFeM8/F2DxMSc3wZ70dDgcHEntlUQowY/bTQCxvTM/A7YqmssD22wuClj/88EP/459//pl+/fpVKf6q0uSulCrXbysyKXS05us33glK7m/NfZjttlbYs+x4HAbnsl/gkrL7e42dWWnDadvlOLL//IXB6U/4t81Kux2HsZMpzUOe22ucZdbt2b2b/z4+jgtvvQPjdYId8pxe0tLSOP6IZpU+n2KHK+T6bbZEnhg7ngEXnxVye3ZeU/Dt63YLm8zBfeO466DLpA4cptQh4Ikpd/PEhCE1drxp44fyxMQ7D/o4hb6EliuxZG3dyqzJdwDg9sYB4PH1EjHl1GJ3Om3k0IINv2WSJS2YmT6MVb8uY9Hnn5BDcrmJHWC7I8r/ePK4DADmPDyTTGnBe/PnkGlPDCqftXYzB2On08YbH3wZcltBwJfCJm/tNqfUFK25K3UI2J+fRKEjofKC5fjsvddZvexbhk+wuvvttrUCT01FB267m7fnPkaOvRkzJwxFSAranmWPq3D/QrEDXmzG8O5rH/m/NMI/v4t5E0ZQ4LS+DDLtZV+r0sn+UGb1Q6ldmtyVOgRUNdmVtvL71ex2JvHHr0s46tgTrJXlNCQ/MfEOotxQHAXGVcSI9PDG+tvvbQZ2DzukFclUHK8xwScvCnh+1X2uFdXy6xuHo/ZTryZ3pWrZ3HH3IbZ8BgW0N1dk+rghOAWGTixJuit//p6vF7zAdltrmpkc7kifBcC88aMolnx2O1sD8OUbL5Ykd58Fz89jy5oscGxnZNosdnqbWQ2ybkBg1tRR3D76UWakDWMXLUnx7mLQxBll4trjLPkpsMMWQ3k/DWalppLjCN10IYbKr142AkuXLqV///61eg5N7krVsix7LBB+O22uvXWZdW+8+THYrPW7JZm5E4YxOH0mmbZEMCXNEZnSnLkTbgexep7MTE21eoU4IMbdNuT5PIXW7fDRXjvYINPWnLnj7sHljQVn6MtyxlY2sW/PzGT2vHlQQa00W3RMKYC9e/fW+jn0gqpq1KaPH8wj42ruQmY43n7tRZ6bHbrf9txx94V1jH3ugC+AUjXhrIAEuiMg0RY6XCEv2u6S1sybMJxMW0k7epY9np3lJPbyzJ43r0rlVe3Smrtq1HJtbbB5q96t7cnpaXh37qX3+edzQv9zyi23/Mcvy6xb/eumctudrVp+ibnj7kFshSClb5+vXtvGTin7qwAgs8zxVX2nyV01et4QTQyVKc4pJseRyPIP3ymT3P/49Vc+eH02xtuWOEc2SEvA+pUwcuLcoMQ+b8Jw8t3JULZLN3Mm3E62PRmIL7Ot9vtaqPpOk7tSFfhh0acs//gTih0eRkwoGc7V48vGO71WTXjWhCHketsTZ88m3msnz94a7B4K3Cn+m19ybW14+9Xn/Md46vEpVo054I7MQBW1T+93upk3/g6wVX7jjmqctM1dqQqs+PQdsuzxiCuBx8em82BqRtD2A138Ct1tcdtduN2tyJKStmtXqeaXzFWr/I+35FU8PGxlMjWx11tSZjjJmqfJXamKGOsj4jF28pyGAoeLx8cNZo9EhSye73BjJHRNHELfUq9UbdBmGaUqVLaG5fK2x+0sqZGnpaWFbDNXqjx1UXPX5K4atJlpQyj2GO6aNLfMttkThoCv98jKZctY+L938doyibU52G88xIuDLJt1MTRfHBwY3rXAXn7NvDKB46UoVZs0uauImjF+CB7j5c5JNdtHevr4wcQ0bc0OWoMdHhs/nLsmltwhOjttCNsDugV++85z5DpbACnkAVEeJ25jwGkNjBU4mqCxVT+5K1VXNLmriNplC93vurRHxg0ngTja9jyGf15+Q9C2JyYM5ohep3HhpdcC8GjqUPY62tBsl4DT6jQYU6q74z5XO3/iBsAbFTQ8d3lDwypVX+gFVXVI+eyj//Ho2EnMTgueijfam0SWPZbMX5YGrZ85fgg7pQ1/LF3PvHF3A2DzzeKzx17y9s5zt2PyuAxmpYaeb7M+jSioVDi05q4OKb8u+pi9zhT2e4LHQfFgx7p1J/hClEOs6vZep4e9+IaBFau2Hnhz0oEuizlYbd7hzBCkVG3RrpCq0fEGvOefGDueJ1LHBW03JpovPnkrYE3wh2RG6gSiTcW30s8ddy8FjtqfCUep8sTExFRe6CBpcleHFAn4f6fTxk6HVTM/UNHJssfx+zc/+MtnlRoTZZdDKu2RkmWP89fulYqEY489ttbPoc0y6pAwffxgbJ5W2GxlB/GalTqWXEdJR/Jse3RdhqZUvaTJXUXEnAnDyPYNqAVQ4D2MQqeLFONgNyVzcwLkOMreIZSWlkayuxi037iqB8445WyM3cXXX39dZ+fU5K7q1EtPPUbe5jUUu9v4uykCeKpxfSlHE7uqY81cNnaXM9BbeewmhgHnnsavv/7qX9ezZ8+aDq0MTe4qLI+OG4RDhJETy97pGWj2hKHs93qJiWlCVFGhf7q2R8YNJdpmsHuSyHG0xmErueMTwGX3PdZp2NQhqoXby/CMNHKzNzN9ztOVlrd5nXhtLq66/CoguIdMixa1P36+JncVlr32FP/j2ROG4DKGkRPnMm/cKAq9iTid29nv8rLf2QancWIKDDsdcaSlpdHUZcPQlp12N3HGAbhx20PfJGRM+NPRKXWwWrq9dD3hKBYv+wuvzWoKbO4yXHDjZbz0ypv+cq08Lk696nwAklofxuHiZqMJM336krrNZvVfOfroo2vwGZRPe8uoKtsurcm1tQGsm39ynbCdVhTbDgOsYW7zA7oa5jm97D9wN2gltfLSMxEpVVqK2el/3MaTX2Z7O2cuKd49YR2rRUoUf7v4WgKnP4m259O5a0/a2/b51w2dlMFxx57qX755wmSauwPezL7d09LSyj1X69bW3djdu3cPK7aDpTV3VaHXnp9LzpotFU56fICrnNp4II+2udQb8S5HyZdyhKV488j1tqTQ4cIEdmN15IOJK1k2cNvY6WRu+ou3n5xNfHJTdmR7ySvVTp7i3YPd4eXft087sFsZt45/hC1b17Fn784QW0vqKUfHeDjjskvYvHplhc+hRYsWjBs3Dru96tM6Vocmd1Whnet+Zocj9PgvT4wdD1WcRLlIx2ypN446vi0/rdwU0RgCa8IPj53sexSQigU6NXGzLs9KZWnpVvmUDkcweNLDgHXTWh4BXwBAQtuWXDu4ZCiK9o5dbPY2PXDIkvXtOgGdKozRGEjp2oeUrn0qfT51ldghzGYZETlPRP4QkbUicn+I7R1E5AsRWS4iv4jIBTUfqqoNU1MzmDF2Qpn1M9OGkjYhDUepgcofHzvR/3hnFRO7ql/+efl/6vR8Do+TNs7ddGpR+fuqDTuwexxcOXgs1985ucKyYrdmvErx7iHJV7dIatkuqMwt4x/nygvPIdnt5u833Fy9J+BTesL1uhhqIJRKa+4iYgdmAecCW4AlIvKOMea3gGKpwOvGmDki0h34AOhYC/GqGlbkcFEUoqnE42oFTihwNw+qApT+easatlZmO9ulVdjlW3uKDuoms8FjrWaSkG3XAW/TwWkzwz7mdXdP5uXHU7n27sls27iObz9+nQsuv7FMue59T6V731NDHCG0JtEudnocNGsWV2pL8OcpUsk9nKrXicBaY8x6Y0wx8CpwcakyBmjie9wU2FZzIaqatOijT5g/fiTTJt0TVvlcnWGo3mph2x1WuWi3E0w1ElDAPs3dVlOJiS4IKmLzOPxlE1wldUnxWqknzh1+y7Dxt8ZULdb4hAQGjptGfEICXY7pyU13ZVS+UxiuvTeV049oxvlDRget79XKmtg8uYn16yCxaXyNnK+qwnll2wGbA5a3ACeVKpMGfCIiw4F44JwaiU7VqEfGDSZRYsi0JdGqKHhy5unjBxMjTgalWxNaaFI/dLQ2OWRLcpX2SXTZGZ4xrcLeGwDxbgf3TB7L1k2beO+paez2NKfA6aZra6vXUnkj8Ng9Dlqxm0y7NRLnOVdeSMvk9rRKSfGfs43JIde0pQhIlmxuz5jN1NQMihwuDm/lJC4+lp0b15BPeM+topSe4t0DtoObcLyqHNExnHXjHWXW/3PYIP4JFBe52PJXFq3b1X6f9lDCSe6hXtPSf/N/A88ZYx4Vkf7AiyLSwxgT9BteRAYCAwE6dOhQnXjVQdhnb8OBzl0uE8W0SfcB1of4QNfGaeOG0ubI2h/USIVvSPos5o0fRaYt0aothxj0zOFx+u8dCErovpvCUswuMqV5uedo16EDgyY+VmZ94IdfvDZay3Z2eVI4tndbtv6cDSSQ6LLT/bgTyuw7OH0WU1ODa8kXXfF3cnN2cNo55wEEjdsf5S15Xm3NLrwENwG2TC4idgdcMLBsQg0Ve6RFRTs58ujDInb+cJL7FiAwwvaUbXa5BTgPwBjzvYjEAC2B7YGFjDHzgfkAffv21WH5atjccfeBbR+D02cFrZ8/I4N92Z6gSZxznYCnbJ/y3fZW7N6YXcuRqooEJmo/36elhduDTTxB4+1Eu52Mnjw2ZC09Oa4Ab1Exg8bPqLQWX5EYt5P7J48NWvdWYSH7Vu0gyplZ7n5xxkURgM1KNT169g3aLgFfH+07lAwnMTB9Rplj3XzHpGpE3niFk9yXAF1E5AhgK3A1cE2pMpuAs4HnRKQbEAPk1GSgqnLWDUCxfP3+R5z+j/P86wt3bGGPM7zp7FTdSXA5uOS6f/HSa6+Fv5OAlNdYYoQm7uAf2rff96D/cYzbSWHprqiVVLGMr4DdlC14yTVDQu7Twu3Bad8FQLcTO7N2+XcMHT87ZNkzr7iNT155i6jYXK4ZlFZxMKpKKk3uxhi3iAwDPsaaZfIZY8wqEZkILDXGvAOMAp4UkTux3i43GRPi3aCq7fFxgynweHHabWC83DN5frllVy/+jOXfLcbp3I7L7SHXkVJuWVU7mrhs7KmgZ1FTyaVd12507tYN8doxvlmjkrxZiGnNrlKtoUW2fCARm30HeJqFPGZaetkurYH+dsEZvPPJZwCkeHeTaQt9nECX3TaWt+c9gTs2/Pbs4ZNLath/u/g6/nbxdeWW7XZML7pl9Ar72Cp8YV2qNsZ8gNW9MXDd+IDHvwGn1GxojdfstKHkFrsYO+VJZj+SRsu2HciztyHaONlf6uf6zAmD8bpbs8tZkgwy7Qm+yZ5b6W1qERItxZR+8cVrw9i8iNfOnROn+9c7vTaKbR6STTa3T5zLglfmsevP4KaOEelz/I9n++aBrewLpLTeJ5/qT+5FUgRAguytcJ+U9u0ZPOnBCsuoQ5PehXKIeSR1ENtpRaxYF5y374Mtq7YAoe/u3CFtghK7qns2r52mrlI1bccemrqCP16xHhtNzA6aR+WFPI7X1wTyr2sG0SKq4qQLYPddWI0yVb/3YET6HJqznXOuLd3CqhoKrdcdYpy+79t8W0liCFU7m5maiq259leMhCRvlr93UVNvFkZsQBv8DdgG7vJdEHx83FBibUKWryvjXekV3XxT8gUxfMyjlV4AjSKftt5CiuOq9zEekRa6HVw1DJrcD1HlDYl7wA6HA/boZY2a1NybTbQ48HriyHZYPYmaejKxO+IZMeExHkjNoNDh4tzLBtO9V3A78azUsYCTJBeIvaS30Z2TZjN3wu0VnjfKGIopuXhZGbHvBlphHHsYmDan0vKqcdLkXseemf0A+ds3YtwtibLtYVBA2+v08YOJtuufpM75+o57JZpB6cG3v985aZ6/WJu2dnIzN5VJ7ABN2jbBkb2VM666laOPOy5oW7Hdic3lIMG+NeTp7VE5pHgNg9IrngjlgCHp1a9xtzDZOPVj3yjoX7mOFW/bbY2y6ABIAmDe+Lsoknz/T/0DZqYNBcIf10MFS3KFvtM22u2k2Ob191Bpw3b2uVI48dzz/WVsXjtNvcHNYTcNLTNmnt/1Q+8rd9uI8dMqjPPOCdPL3Sbemh1FcHi61vQbC03udSyXRAKnl0tLSwNbE8SbBOIJKrtDE3u1paWl8e3Cd1n55Vf+W+QPGO27GSd9/CSMzVPmpi+A8RPH1UmcFWnuzcZL3Q0RqxoWTe51aHrqEIrKGRv9QC1SlWhttpMdxoiEzd2GXY6yPYZOOfsiTjn7IuaPH0GhFGHcKey1lVTlE002xnPovu4jJmotW1WfJvc6lGdrBxwaM9vUB0PSZ1faY0S8NkZMHs/DYyf7Zw0q3QVx4MSyt7ID3DUp/GFjlapvNLnXkDkTbidbkklw2YlxZoK7BTscdtp48tln8tjnSNG7Cg5CkltwyTbyOQyv3Uri4rXRu9eRAETjYj9CijePQRmPRzJUpQ4Jmm5qwPxZU/xDsu5zethBK3Y4rLbSLHucldhVGc1cZZtSHB4nKd7dxLiDr4TGtDTcPWkesWYrLUwWbTz5JJHJRZdZt7a7JJumLiFfm7eUArTmXiNM9nYIY5yOxqqVu5jtDmvEvxSzkyLjBvFicwjQhpbebFyeFP8sT4Mmlu1dMuiONADumRy6HXrU5PC6ESrVWGhyP0jzxt8Z1gBMDV0rs50Cdwp7nWVrznb7XsCasODAZCClzZgwGOsuT6VUTdBmmYOUaWsa6RDqjhFaustv9nDIlpDrjz/zcpLdLlqarHL3dfkmaojx6hytStUErblX04zxQ9hlazxjpNs9DqK8W0FallPC4C1n/LITB5zBiQPOqPD4oybNY+6EYezV3kRK1QituVdDxr23EWUSKi/YgMSazdyXUX6/a6/AyWdeTQu3l5be6s3kNDh9JvdkaNu5UjVBa+7V4IprR/kNDA3T3ZPnlbvN6XbS98wrOOnsAZx09gCAg5rSTSl18LTmrqrIahtv49nvXxNnNtHvnLOCSsUXZ1kDch2kFO9uWhudsVGpqtKau/Jrw3bwxJNljweghctLtH0v+02Bv0yBdwfNXa3ZJ/mkeAw5JHFnRtla/T1T5vJExn14XAfXhh6qW6RSqnISqalO+/bta5YuXRqRcx+MwNvc65MUs5NMaVGy7M0li5b+MW3S7hkO0YngiPI3qbRweRmeMTES4SqlyiEiy4wxfSsrp80yVTBl9MB6mdjB6l/eJqB5Y9DE6Yiv2aSpKxPiW4DvRqNWxrogGmXfXfeBKqVqhCYlAgLdAAAgAElEQVT3Kmhir39D8No9DlK8uQAMTp9FvMtBc7fv15qvSdxTqml8aPoc4j1ZDCpnwC2l1KFP29zDNGfC7exwJEc6jJDiXA7yQ/yiSPHmMmhS8EQQ92SkBixZSV5sZZvm7pmkXRKVqs+05h6mQvchfMOSbAtajHVtJ969LWgKv1ASPVZSd3n0baBUQ6M193qsqcuGzb6NkZPn8sLT01i/2Wojvy8jvDk2XWTT3OtlhN44pFSDo8m9nmnmEjzYiHJsZXhAUr7hljuqfOPQvRllp5dTSjUMmtzDdGA42rrU0mSxQ6yREpNNNg4TR0yHdtw46K7QOxhI9rhCb1NKNSqa3COgqctW6ZdFgiuTYRnzmDF+CA4RhoYxa31aeloNRaiUqu80uYfhhWdqdqJim30blY1dfrfvrk+dJFkpVR2a3CsxffIYct1RNXa8A+3iMx+egHPfLow3lix7HHaPA49vbtCaGJNFKdW4aR+4SsQU76uR4yS7XaR49/qXh92TzqD0J/DarXWJ3pJJMBLN1ho5p1Kq8dKaeyWE8GvRTV02nFLMsMmTAWu4guLotjg8Dm6fnBZyn6Fpc5g3/g5yvYUcaKoZNXH+QUatlGrswqq5i8h5IvKHiKwVkfvLKXOliPwmIqtE5JWaDTNyTBgvUZJbaJscxZ0Z4/2JHWDM1Pk0se2jSWJhhfsPmjiN+6fOBWOjqUt/TCmlDl6lNXcRsQOzgHOBLcASEXnHGPNbQJkuwGjgFGNMrojUv0FYylVxso3K38rIh54sd/td4x8J+0xp6ePDLquUUhUJp5p4IrDWGLPeGFMMvApcXKrMbcAsY0wugDFme82GGTn7Pc3K3eZ0OxlTQWJXSqlICSe5twM2Byxv8a0L1BXoKiLfishiETkv1IFEZKCILBWRpTk59WN2nT1OT8j10W4nUbI55DallIq0cJJ7qCuKpYcRdABdgDOBfwNPiUiZKq8xZr4xpq8xpm9y8qE5wmKgqaMHl7utqezTkROVUoescJL7FuCwgOX2wLYQZd42xriMMX8Bf2Al+3rN7Whf7ja7VHyRVCmlIimc5L4E6CIiR4hIFHA18E6pMm8BAwBEpCVWM836mgw0Evw3FZXS1CUccezpdRyNUkqFr9LkboxxA8OAj4HVwOvGmFUiMlFE/ukr9jGwU0R+A74A7jHG7KytoCPNbs/ib1dcGukwlFKqXGHdxGSM+QD4oNS68QGPDXCX71+DluTNYsREbWtXSh3a9A7VKkh02RmpE1sopeoBvR2yCkZljIt0CEopFRZN7kop1QBpci/HnAlDIx2CUkpVmyb3cmSXGh6ntakfd9QqpRRocg/pwTG3lVk3JF0nk1ZK1R+a3EMoiCo9dI5SStUvmtyVUqoB0uQehpbe7EiHoJRSVaLJPQzDJs6JdAhKKVUlmtyVUqoB0uSulFINkCb3yhh9iZRS9Y9mrkrYvfoSKaXqH81cpUybNDrSISil1EHT5F5KrKsgaLmV7IpQJEopVX2a3EvxeOOClne49kQoEqWUqj5N7qVsd0T5H9s8DsZOfSqC0SilVPVocq9AlJFIh6CUUtWiyb0CscYV6RCUUqpaNLlXINqWH+kQlFKqWjS5B3h4/Kig5biOR0QoEqWUOjia3AMkGnvQ8g23DIlQJEopdXA0uQfIssdVXkgppeoBTe4+LzwdPI2e6LADSql6TDOYz95Nq4KWkyUrQpEopdTB0+Tus8vbNmh5aPrsCEWilFIHT5O7j8fujnQISilVYzS5K6VUA6TJPYQWLm+kQ1BKqYOiyR14/umZQcvGnhOhSJRSqmaEldxF5DwR+UNE1orI/RWUu1xEjIj0rbkQa1/WX3lBywVubX9XStVvlSZ3EbEDs4Dzge7Av0Wke4hyicAI4IeaDrI2TR09kAJH8ABh9015MkLRKKVUzQin5n4isNYYs94YUwy8Clwcotwk4CGgsAbjq3VF0W0rL6SUUvVMOMm9HbA5YHmLb52fiBwPHGaMea+iA4nIQBFZKiJLc3Ii365duq1dKaUainCSe6gZK4x/o4gNeBwYFaJc8E7GzDfG9DXG9E1OTg4/ylpSuHF9pENQSqlaEU5y3wIcFrDcHtgWsJwI9AC+FJENQD/gnfpwUdXrdZZZ18LtiUAkSilVs8JJ7kuALiJyhIhEAVcD7xzYaIzJM8a0NMZ0NMZ0BBYD/zTGLK2ViGvQTmlSZt3wyZMiEIlSStWsSpO7McYNDAM+BlYDrxtjVonIRBH5Z20HWJvcdp1GTynVMDnCKWSM+QD4oNS68eWUPfPgw6p988bfAbZmkQ5DKaVqRaO9QzUzRGJvqTcvKaUaiEab3EOxNw/rh4xSSh3yGmVyf3T8wDLrmrhsDLkrre6DUYeMlcuWMeWRf/P91x/V6nkmTf4PUx75d62eQ6lGmdz32srelRpvz41AJPXX5NTBkQ4hpIwHr+ehGRdVa993v57KSb1/ZPGqR0NuX7zoM56YfyYZaQc3cfrJJ3/FSb1/PKhjHIwpD19DxqRBVd7v2y8+ZskP39RIDN9++T6PzLyAzC0bauR4gXJ3bWfXru1B6/bk7uTJuWNr/Fw/b1jCJx914Y1Fh97kPo2uHWLq/QMhpmxy79D91AhEc/AemnERXmPn/pFvVVr2medOYH9xPMMHfnlQ58yY/B/6n/UVD8/8B6d1u5vv/5zCvsyOjJsYuTF5fl78A+8vmka/E6ykufDzTqxcczgjB30e9jGio/cCkBCzN+T2zxc/xUm9N7Mtdk+ZbZ9+0pnf/upQpfOFkpF+K/1O+4Ki/Q42ZCcz9NZFB3W8MsefNIh+p/yA8cIjs87HYXdxx+DPyi0/dfJgvM69jL3vZQrNUHZlOoHfKz3PQ9MvplXSVm66IXSP6J/XPMDx3bfx0pvDuWf4u6xb+zuZG//ix58/xxm/lChHIYNu+b5sPA/fQKfDltPz2Jc46pieIY/99RdnkpBUxLYtt9G3z9kc1uFoPv30DI7sWsCsmYUUeHJJabmck098kSO69Ajad+HnH+J0ODn99HPKHPe1V49DgCuv/sW/7rulc+na3Ev+xmfg1KFB5d3uItxF+4mJb47Hk49IFDZb3aXcRpfci0IkdoDzr/5XrZxv8rjbOO74Rfz6W3/Gpj7jX//wExdStL8Zqfe/VPH+YwcidsPYchJnnx6/hR3L4R12Abt4dNZ5FGV2YszkWWXKTEy/DbvDzdixz5Z7HHvCbgBSmm1j4dIn6XfCejITd/q3PzvvUbZs+5P27Y+ibfs5/LzmaO4d8S4AS778mh37bqW42IGI4dc/+2Db14TRk2bz1uvPsnHDn/To+zrLV3fl7ts/DPu5ffxjBif1Dp4Ht0eXjUyecjOIl9TRz/vXT06/lbioJtw1+jFmPXk6e/e14v4738D47rsufUv2c/MfZ+u2NUiiNc5/04QCJj9wI2JgrO+4NoehR5eNZeKa8si/seU3x20rpGnyev5+asmQF2+9/iyXXHlzUPkTT/4CgOh4N0cdmek/f+cux3LqgLIJ54BHZp3PcZ3W8uPSM+jc6UfWbjgeccfQ/vClGCP07zuPo47pQ98TrC8fscHx3f4sifPRqzBuJ2Pve4ln5z5E0+atuezKGznx5E8ByHjkWvr1hrgmJd2HF37eid/XteP2274G4L13X+bPzc/z99MfoM+xKwFY8H/HkJd7DUUF+SS2+ICUtnvo2eN7YpzWEFS9j/mNtxZ0JzGpCAR69ip5Tq++0pOcfUnEOAuxi5cLz3uLHl1+JK6Ji/c+epJVq3uydtNPdD3iRC651Hod9+7ZQ0JSEQBt2z/Jtuwn2ZYNzVpax9zn3kzf45YBsG7jxfyxxo4zxsOK37visB1Hj65v4HLDI7NP5fijrS/W43ouJrlFMi1b7bf+bu9MZ+v2d7n91rJfil8vns1h7U+kaUIblv0wAJvTy9lnrePLr46ledIpHH/8C+X+DWuaGGMqL1UL+vbta5Yurfv7nNLS0sqsS/HsY9CkR6p8rMkPXs+JnS/nb/8KNY6a5eGZ/6B399/ZlxvNr2v7kHrfi4D1wQD44acT6NzmHHYWPsv+wqYc03EtS78bwLgp85j5WBrdelnl9+yK5vfN3XHtaoY9oZCx970UdJyzz1pHeupAnEm7SU7aRM76Ppx01ke4i2z8/fw1QWUBtm1tyvXX/8T8Z/rjtLvJXHMiACedZbU3L/7+dFql/E60s4gNf/Ymse1GunfcwO8bO9CpzVZim7jYnxfFyvXHcdLxS9m5PZ41v50E8fs5oeeP2KMM339zNv1PW4jXLZz7t7WAldz3eIMTWnZmIm6PnXbtd7Pk516c0HMFrgI75/3jz6C4v/vmLMRtJ3XSXACmjB9CXNs1HNv1r0r/Vt99cxbjJjzpP97uHbH868qV/mOv+as1XY7IBsDjFuyOsp+LH5b35aTjK37PLvm5J+49zbn8spH8ufZK4psWB21f/M0A+p32RdC6ZSu7U7y7DVHNsunTI/gLqmO7D9iw9QIAjNdKygA/LDqXMePn8uC0S3DvbclJ/b/wbwslKzOR3P1N6dZ5S5ltW7Yk0b691Sx59lnr/K/J9z+eQv8Tvw15vF9+78RxR6+zYlnelxZRR+NM/JDDO+zkp1Xd6X1M+ZWO3Jw4kpLzyw+2HDnZCSS33gfAug2t6NQxuOll6fJzSGzyG0d12hZq92r7+c9z6NzuWOLjH6+wnCfpQey595VZnxw/mJz91nv2mO6P07r1P7AG260eEVlmjKl0BIBGldwn33cr7tj2ZdY77JvB7iV1TEltdeL4Wzmi4/G0SmzLr9lPsX9LChMeeIrJ6bdiDCS02uR/c3vdwg/L+2NzFuB0FtH7mN/47qtziGu9kcOab6VFq5I38uKl/WnX+ncOO6z8Nv5Vfx5OfmESrZpt8tW2Q/v+6wH0P91KFIt/6ke/3ovLLbt1azPatdtd/otTTa4CO85YD+5iG8YrOGNCD9+wem17XLv60b3n/3BER2amq+9/7E/q/S8Ffckd6rIym9AmpWwzEEDOplEkdwh9faC6flt7GN07b668YCOy/PeTaNNkNSltQ/8dquP005bjdJa9Qz4cmtxDSJuQFnIYtNNOt2rHZ59lJev33zmKmASrz3vezhiatihk/YZW7M8awLH9XjuoGDZual5hwlZKNXzNmp5Anz6vVmvfcJN74+otE5DYTzv9RXoc+ylHHbnQv+7hmf8A8Cd2gKYtrLbBIztuP+jEDmhiV0qxO29JrZ+j0V1Qbd5iM82aZgGQlJQFSSXbenf/nQ/f60pUXISCU0qpGtKokntcXC7HHPNlhWWi4nTIX6VU/ddommVmpw2lT98KJ4pSQBZtuI/H2EtipENRSh2ERpPcC6t3YbrReZdL2CKHs4STIh2KUvWGFxtzGM5fHBHpUPwaTXLv2evjSIcQ5F0u4RdC32EH8D0ns5Lj6jCi0kLNrqjqmgHWE17XzSKieIdL8TSej3WVeLFR1b6Bv9ONnbSotNwgnmWRnMl07qmw3P+4gocZXcUoqkffBTWskOhKP1xebLwq1/OghBwSH4CZMoqpMuGg4/mKATxF+OPAhErpsxnBLbzELpqH3OcH+vM73aoVnxtHua+XGwfuenBZqIAYskgJWvctp/EdJUNa7COBIqKqdNxdNOc6WcA4eYhvOZXttMaNgyzahCy/gKt4Ta5jEWdU/UnUkj84io85v9r7b6ID6UwOeu120pyBPMciTmc7rcKqLbtxcL38H69xbZXOP0kmM4onKi2XLwkASMDXRy7NcJV6/y6Qq1lRRzOQanIPg5Vkgu8oW8zJXCsL2E6rgHJ2bpFXmMEo/uAo8on1bzPAl5xFMU5e55oaietdLmEVPSosM1+G8YWcyzo6UUwU+cTxLadxL9MoJBqA7bTiff6JCwe7fd2HAms438oZFEosw+XJoPUrOZbVdGeG3M0kmVxhHLtpxrWygE/5e9D6G+U1xvAIxnfOA0lyJcdyo7zGEJ4pe7AatIX2fM/J5X7BFBLNXIbxNpfyI/3Ioyn4Yr2D2XzG35jKBEbJTAqJ4VpZwJtczmy5g1lyJwDFOBkkzzOGkhuOskjhF3qFOiXFRHGtLGC4lAw5MVvu5E6Zzf08xiiZRR5NyvzyK/C933Jo7V+3lwTe5lL/3y2PJrzJ5UF/xx/px71Mw4uNQqJ5hesp9iXT/cSV+eICGMLT3MCrXCsLWMKJuLGTQ9lJ7yfKFF6QWwHwIizlRFbQmw+4CAM8y61cKwtYTfeQX+QvcxN/Sjf+5Gj/ui84l/2SyBwZyZ0yh1R5hE10AKwvxJe4iWtlAZ8FvNcOPJ935TK8AX/rfSTgLedX6oHPtkuiMRDWL6IDyd0Aw+RpZnFnyHLHHTe/0mMdrEO/WlQDHp5wM70DKjPFRPEa13IJbxBFEf+R/9LJ/MlERpNFCq9zDTfxJB7seBFG+D5kLxtr/Jm9JPrfON9yGpm041spOcFS6cdS+nGEWctgZrKTFrhw8qTczi+mF0tLtWfPZRh9WMIajqIVWSwMeFOuoQsuouiOdVu6CweTmchFvMWrcj0AQ800ZssdTDT30ol15NKMrXQIqvmPl4c4xXzF95yK13fr83vmUt6UK/xllpiTWCPWh8iU84afwwiGMgOAqZIWtM2L8BedmMAUjO8cp5vPuYGnuV2eBuA5GYjTuMmmNV9gjZWyRQ7nOhbQw6xgpfRinEn1Hztf4nnC3MliOZVO5k/uZxJx5LOPBAbJ88SbvczlZv8PbjcO9hNPU/Ks/YllLV1ZSU/el4tJNtlMwxrg6RVu4H2xho64zLzGMk5goxwJwDiTypGs5WMu4BsZEPQ8W5gcdoqVyJ6lZHTFW+RlAN6QkuF83+UStmLdFZ0lbf3fmqOkZIyZazZ/xNntX+Qt/sUqjmV7OTVzgExpB8BQse6m/od5i884jwzu5ltOB+BNuYI3uYJeZikO3CyVfrzOdQD0MktZIX3pZlaxns7sI5G3xXpfv2cuxoWD9+USXMbJjTzDBB4gU9rRyyyjB7/Qn0XsI4E9UjLZzTMMIp943OLkNjOLjqynKbv546dO0Mcqs4oe/E53/idX+ffbaDqySM4EYLJMIsHsYS438zj3cS4f8jSDyZGSL6pckniIVIp8lZJAo+VxXjb/CvpCfFYGcrb5mB/pxzyG+dc/tO9xuvzlpr3jbWZ0H0ms2U9rstng+9vHmHwKJY4WJse/z3WyAICLzQKiKeR1uZYLzZsspy/9KRncLVtSeN7cwlKs4TyWSD+KTRRu7MRR4C/naFr7v64axR2qb7x6PEmtSm4dfo+L+a/cUCfnrsytZjZPydBKy40wjxBNIQ9Lah1EBd3MSlZLD38yqG+amVxO2LCPT484LNKhKFXG0MNaMb5z6EEMKxPuHaqNouZ+ILHvJYHB8nwlpetWOIkdYIbcXcuRBFstVnNPfUzsALsliU+PSKq8oFIR4K3ypd2qa1Rt7odaYldKNU4FntofPK/RJPevGFB5IaWUqgPPb9tZeaGD1GiS+3wZVnkhpZRqIBpNcldKqcakwSf3D1/5v0iHoJRSda7BJ/dlvyws9yYFpZSKhBhb7eekBp/cEzquZT2dIx2GUkr5fdina62fo8En92O7/sUEeSDSYagGotPuokiHoA4xx23fX6XyA9sn0y0htvKCB6nBJ/dwnbsmG6en7I0FV/20JuxjvLhrEyMWvhNWWbvX0KT44Pq6Hr2roPJCPjct/iVouXlheJOSDPr6O65YsT5o3fBPX2XYwvcr3O+5fWuxBdz9fMuipWQN6MWIT14KM+Lw9cnay/U/rGLIV1+XW6ZX9j4uXL2VPpl7Sf/lg6BtQ79YyBkbQndNG/bFR/7HVy5fyy1bv+blgoonkD53bXaF20d89malr//Z67aH/V46789MRix8i0FfLeLcNeWfu+Oe0F9M5/+RCcCwzz/hpu9/DtoW4zbc8MNK/3LWgF7c8fm79M0sO1n02eu202NHfpn1B7Tf5/I/PnGbNTxE4Gdu4DeLGfLVN/7lq5av5YoV6zjvzyyOz97HoK+/B+CWb5dx6cpNdNxTDMAJ2/Zwx6Ivyj1vdfUMSNpXB3wGTtxmPffeWXvpk7WX9y87kY+TEpi2cx9ZA3r5P9crOrUha0AvVvftxqk4/ftPqOadqVXV4IcfWPh5J671jQsRaPinr/J7xwF0W/8pYx54xL9+ypixzDi3ZLyVrAG9yEjN4M2T/saWBCdXrFjP//U6MuhYTo/BZReyBgQPBDXpvpF4Y47glF6dOefSC8lITWd3Qmcevr9kZLqp949m3WGn0N/2O2OPPsd/zuGP/4+4wv3EF25n9pln0zrfzc//6Mtts97j3e7WWCXDF36I2xaN07sH483HFXUYTnceuxOOYF98Iss7tGFDE2vApJEfP0d+7NE8eXo/bvr+Zx4YcyNDZrzNm8ceTsc9RVzyw6t4bEm4bQnkNmnPtLuvLPOaPTj2Lly2NqROuheA9LQZzDnj9KAyp23axf/deJb/tVzX/hQuaJrNv6692V+mzRcrOHVTLr3/eA8kgf+echE5sdbN0kO/+BSHezc5SccR7Srkuf7W4FjDP3udsRlT/MeY+XA6k/ta48KMXPgKoyc/5D92oKt+WkPrnSsYMyUjaP3ECdPIj06iyf4/GZORwaplPzFu1R7iiorZ2KIJa5JiGPrlQsanj2JMxrMkFf7OPZMeDHoOYH3B/t48llsXLcFrcxBfsJWxk1MZO+kpYly7iLXvpsDTjD8O68VnnVtxxYp1PHHnv5j+YBpTT7yEmxb/wnP9rKGdj8wr4qrVH7BVuvDQ/dZYMBnjHmRLq260yVnLpjZHM8Cxnm+LOvBmjw4ct30/fdauYerYm4Ke26b16/lj5S+sXP4Tez1JuB2x5Ca2on/8dh5J6cO2BCejvvwfHx/zd47duIXH77kqaP93X3+Rpav30rdbIj169eaIrsfQ5osVXPrrRuaMsF7zlStWsPD9Dxg5dkyZ98kB7Rcuxx3Qtjzyy/dJTo7n89hjOHLbBp469QR6Z+1lkH0jP/5ZREbqLUGvbenPE8D0KY8ycsyokOc7sF+bfDdZcdb7afCiH2gW4+bbVt2ILyzio64p3PL9ck7r1YrvVmQzv39vAC74fRuxxcXEFRWwuWUrmu3fz43dEhlYnEhOrIOsAb2497G32BvlYM6wC1m9+i+6dQs9GuVfqzfwx19ZnHdBP/+673L3cdmKtfwjuSlP9zi4Md/DHX4grOQuIucB0wE78JQx5oFS2+8CbgXcQA7wH2PMxoqOGYnk3jrfTbzLw/qm0Yxc+DKjJz9c7n6pk54ktngXYyfdB8ADo+/EY2/G2MnWYFwZY0dj84DXkUBifDF5u3cz7sHpBxVrqDf1lNQJzDj7Utrud/HThSeQMXY03/a4gFN/fjvoS6k8ozOe44NePbhn2w9cd9vtQdsyxtzBppSz6bLzG+5Oe6haMd8/5XliiveQm9iG13p34ebvVzB1zE0V7pMx5l6aNo9j2N1p/nWzH5lI3q7djJ7yWFDZl5+dxcY128ok5wPHibLDPZPKxv7ivCfYu28PQ0eNrfJzunX2+7zXrR0jP13A6CmTQpa55unP+PzIlvzn25945pTejPzsv4zOeDBk2co8MG4MbhJJnRTeON/3T3me5/r35D/f/sSU1P9U6VxPz3yUXTl53JM+sTqhVskrT81hoetwMPB+t7Y8+Ns33Hj7cADmTnuUtJ5nM/Dr75k4YUjQfhUl94osePlV3C43ToewdsNOrrzqAjp2KbnetmfPbv5Y8zsn9LGS7v9efpmhbY9h2Hc/kjp2YMhjrvhuBevXbeWy6/9RpVhKM8Ywe3MO16Q0J8l5cKO+1FhyFxE78CdwLrAFWAL82xjzW0CZAcAPxph8ERkCnGmMuSrkAX0ikdyvXL6O1juWkNekh79mdCgZ9MQ7tMjbEfSBzRhzP0+cezUnb97N/244M3LBVeLHz7/ko6+Wcd4ZfTjxrDMjHc5BefW52axbu4uxk8sfpO3zD9/j12U/MjK19pNkaRvWrOGFF9/jhusvpGOXLnV+/qrasnE9Sxd9yyXXXh9W+eom98aiJpN7fyDNGPN33/JoAGPM1HLKHw/MNMacUtFxI5HcR3y2gDEZoWtih7L0tBnYitYxburB/TJQqj6YOuEhEmJh+P33RjqUQ1JNjgrZDgi8grQFKpxg8xbgw3KCGggMBOjQoUMYp65ZXk/VrmofKiakjYh0CErVmdHpmtRrQji9ZUL1tg9Z3ReR64C+QMjGbGPMfGNMX2NM3+TksrO21LbUBx6rvJBSSjUA4dTctwCBMx60B7aVLiQi5wBjgTOMMdoZWCmlIiicmvsSoIuIHCEiUcDVQFAHXF87+zzgn8aY7TUfZvWk339rpENQSqmIqDS5G2PcwDDgY2A18LoxZpWITBSRf/qKPQwkAP8nIitEJLy7L2qZPdYWctJdpZRq6MLKfMaYD4APSq0bH/D4nBqOq8bsJz7SISilVJ1r0MMPSGwee2kS6TCUUqrONejk3q3jb6yiR6TDUEqpOtegk3vTlvlsxBrH4cLVWyMcjVJK1Z0GndwBvpKzAfDYGvxTVUopvwad8QoLSy6muhz2CEailFJ1q0En95iYkuEGOm35M4KRKKVU3WrQyT2Q050X6RCUUqrONJrknjq56uN6K6VUfdVokrtSSjUmmtyVUqoBavADrySaPBzFCZEOQyml6lSDrbmn3WuNCLlXmpIbrd0glVKNS4NN7jFtcyIdglJKRUyDTe69u/1WeSGllGqgGmxytzu9eH1P77JfNkQ2GKWUqmMNNrkDbKMtAHnxOqa7UqpxadDJfSppACzsVPeTcSulVCQ16OTuwAWAGBPhSJRSqm416OS+h6YAHJ+9L8KRKKVU3WqQyX3M0JsBiKIYgGP/2hDBaJRSqjDtG3oAAAhqSURBVO41yOTea8B6APZJIgDxhVmRDEcppepcg0zuLVpsASDWWOO59z36yEiGo5RSda5BJvcDCsTqAnnBNf+KcCRKKVW3GlxyT0tLi3QISikVcQ0quU8fm8bJp/w30mEopVTENZghfx8fN4jO/X/EbncD8BvHRDgipZSKnHqd3Gc8mEpU/m52elJo1zmTuLg9/m1PcFcEI1NKqciq18l9V4Gdlq32keDewGEdVgVt2yPNIhSVUkpFXr1N7g+kZnDaWS+F3BY42MDwzz+GAb3qJiillDpE1NsLqiec9Uy520a6X/Q/HjvpvroIRymlDin1MrnfO/jykOvd2LlWFrDTGVfHESml1KElrGYZETkPmA7YgaeMMQ+U2h4NvAD0AXYCVxljNtRsqCX+fuVyALbRjntkRrnlLv/5L22SUUo1SpUmdxGxA7OAc4EtwBIReccYEziP3S1ArjGms4hcDTwIXFUbAQP8jytYIFdXWq5dziLg0toKQymlDlnh1NxPBNYaY9YDiMirwMVAYHK/GHwzY8AbwEwREWNqfiD1QS/M5u3DKk/sWQN6aa1dKdVohdPm3g7YHLC8xbcuZBljjBvIA1qUPpCIDBSRpSKyNCcnp1oBt84LHuHxqF2FAFy1fC0jP36NjDXfWIldKaUasXBq7hJiXekaeThlMMbMB+YD9O3bt1q1+onDJzIx1IYBvYDQF1qVUqqxCafmvgU4LGC5PbCtvDIi4gCaArtqIkCllFJVF05yXwJ0EZEjRCQKuBp4p1SZd4AbfY8vBz6vjfZ29f/tnVuoVFUYx39/tIRKU9NCorwEBfqSJx+M0oeC1NPFKIqiB8kgughZ9GAIEr2pBBFFUiSpWJqV4IOSEmEPpaJ2vGV6jpcgOx3TCnvoZn09rHVyn/HsOTNz9p49M3w/2Mzim7XX+u9vrf3NXmvNzHIcx6mMAadlzOy8pAXAp4SvQq40s0OSXgF2m9km4F1gjaQuwhP7wCuejuM4Tm5U9D13M9sMbC6xLUmk/wAeylaa4ziOUytN+QtVx3Ecpzwe3B3HcVoQD+6O4zgtiAd3x3GcFkRFfWNR0k/AdzWePgY4k6GcrHBd1eG6qqdRtbmu6hiMrvFmNnagTIUF98EgabeZTStaRymuqzpcV/U0qjbXVR310OXTMo7jOC2IB3fHcZwWpFmD+9tFC0jBdVWH66qeRtXmuqojd11NOefuOI7jlKdZn9wdx3GcMnhwdxzHaUGaLrhLmi3piKQuSYtyrus6SZ9LOizpkKTnov1lSackdcSjPXHOS1HbEUmz8tQt6aSkA1HD7mgbLWmbpM74OiraJen1WP9+SW2JcubF/J2S5qXVV6GmmxJ+6ZB0TtLCInwmaaWk05IOJmyZ+UfSLdH/XfHc/jatqVTXcknfxro3ShoZ7RMk/Z7w24qB6k+7xhp1ZdZuCn8bvjPqWq/wF+K16lqf0HRSUkcB/kqLD4X3MQDMrGkOwl8OHwMmAZcC+4DJOdY3DmiL6eHAUWAyYb/YF/vJPzlqGgZMjFqH5KUbOAmMKbEtAxbF9CJgaUy3A1sIu2ZNB3ZG+2jgeHwdFdOjMmyvH4HxRfgMmAm0AQfz8A+wC7g1nrMFmDMIXXcBQ2N6aULXhGS+knL6rT/tGmvUlVm7AR8Cj8T0CuDpWnWVvP8qsKQAf6XFh8L7mJk13ZP7/5t1m9lfQO9m3blgZt1mtjemfwMOc/H+sUnmAuvM7E8zOwF0Rc311D0XWBXTq4D7E/bVFtgBjJQ0DpgFbDOzn83sF2AbMDsjLXcCx8ys3C+Rc/OZmX3BxTuCZeKf+N4IM/vKwl24OlFW1brMbKuF/YcBdhB2PEtlgPrTrrFqXWWoqt3iE+cdwEdZ6orlPgx8UK6MnPyVFh8K72PQfNMylWzWnQuSJgBTgZ3RtCAOrVYmhnFp+vLSbcBWSXskPRlt15hZN4TOB1xdkDYIm7Ykb7pG8FlW/rk2prPWBzCf8JTWy0RJX0vaLmlGQm9a/WnXWCtZtNtVwK+JD7Cs/DUD6DGzzoSt7v4qiQ8N0ceaLbhXtBF35pVKVwAfAwvN7BzwFnADcDPQTRgWltOXl+7bzKwNmAM8K2lmmbx11RbnU+8DNkRTo/gsjWp15OW3xcB5YG00dQPXm9lU4AXgfUkj8qq/H7Jqt7z0PkrfB4i6+6uf+JCaNUVDLj5rtuBeyWbdmSLpEkLDrTWzTwDMrMfM/jGzf4F3CEPRcvpy0W1mP8TX08DGqKMnDud6h6Kni9BG+MDZa2Y9UWND+Izs/PM9fadOBq0vLqTdAzwWh+HEaY+zMb2HMJ994wD1p11j1WTYbmcI0xBDS+w1E8t6AFif0FtXf/UXH8qUV98+VunkfCMchG0BjxMWcHoXa6bkWJ8I81yvldjHJdLPE+YeAabQd5HpOGGBKXPdwOXA8ET6S8Jc+XL6LuYsi+m76buYs8suLOacICzkjIrp0Rn4bh3weNE+o2SBLUv/EDaPn86Fxa72QeiaDXwDjC3JNxYYEtOTgFMD1Z92jTXqyqzdCKO45ILqM7XqSvhse1H+Ij0+NEYfG+xNXO+DsOJ8lPCJvDjnum4nDIP2Ax3xaAfWAAeifVPJDbA4ajtCYmU7a92x4+6Lx6HeMglzm58BnfG1t5MIeDPWfwCYlihrPmFBrItEQB6EtsuAs8CVCVvdfUYYrncDfxOegp7I0j/ANOBgPOcN4i++a9TVRZh37e1nK2LeB2P77gP2AvcOVH/aNdaoK7N2i312V7zWDcCwWnVF+3vAUyV56+mvtPhQeB8zM//7AcdxnFak2ebcHcdxnArw4O44jtOCeHB3HMdpQTy4O47jtCAe3B3HcVoQD+6O4zgtiAd3x3GcFuQ/16B1oHL2XwMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "continue_train = True  # set to True to continue to train\n",
    "num_big_steps = 50     # number of small steps\n",
    "num_small_steps = 200  # number of big steps\n",
    "if continue_train:\n",
    "    for k in range(num_big_steps):\n",
    "        for j in range(num_small_steps):\n",
    "            dXY_list = np.append(dXY_list, np.zeros((rep, 1)), axis=1)\n",
    "            dX_list = np.append(dX_list, np.zeros((rep, 1)), axis=1)\n",
    "            dY_list = np.append(dY_list, np.zeros((rep, 1)), axis=1)\n",
    "            mi_list = np.append(mi_list, np.zeros((rep, 1)), axis=1)\n",
    "            for i in range(rep):\n",
    "                minee_mine_list[i].step_minee()\n",
    "                dXY_list[i, -1], dX_list[i, -1], dY_list[i, -1] = minee_mine_list[i].forward_minee()\n",
    "                mi_list[i, -1] = (dXY_list[i, -1]-dX_list[i, -1]-dY_list[i, -1]).copy()\n",
    "        # To show intermediate works\n",
    "        for i in range(rep):\n",
    "            plt.plot(dXY_list[i, :],label='dXY')\n",
    "            plt.plot(dX_list[i, :],label='dX')\n",
    "            plt.plot(dY_list[i, :],label='dY')\n",
    "            plt.title('Plots of divergence estimates')\n",
    "        display.clear_output(wait=True)\n",
    "        display.display(plt.gcf())\n",
    "    display.clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuously train the model with MINE. The following can be executed repeatedly and after loading previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXecFEX2wL9vlxwlKZIkCCoqgq4RTGBAVMwnenpmVEQUE1GQjHImFEX0MCtiuDv0hwEV9cygogiiIpIkSBBB8u6+3x/duzs7O6Fnpmdmd+Z9P5/9bHd1hdfVPa+rXlW9ElXFMAzDyA5y0i2AYRiGkTpM6RuGYWQRpvQNwzCyCFP6hmEYWYQpfcMwjCzClL5hGEYWYUo/ixCRD0Tk6jTLsJ+IfCMiW0Skn4f4d4nIc+5xCxH5S0Ryky9p5iIifxeRd9Ith5EeTOlnGCKyVES2u8pxrYg8KSK1YsyjpYioiFRKgoh3AB+oam1VnRhLQlVdrqq1VLUgCXJlJKGepao+r6qnJKm8tDcsjMiY0s9MzlTVWsChwOHA0DTLE8g+wIJ0CxFIkj5uhlEuMaWfwajqb8CbwEHB10QkR0SGisgyEfldRJ4Rkbru5Y/c/5vcHsPRIrKviHwoIn+KyHoReSlcuSLSU0QWiMgmt+V3gBv+PnAi8LCbb7sQaVu55WwRkVlAw4Brxa1WEeklInOD0vYXkRnucVUR+aeILHd7PJNFpLp77QQRWSkiA0RkDfCkG36HiKwWkVUicrVb1r4x5HerW5erReSKALmqi8i9bl3/KSIfB6Q9SkQ+devqWxE5IUK9NhGRV0VknYj8GmgeE5EjRGSuiGx25bsvwrO8XEQ+DkirItJHRH52632UiLQRkc/c/KaLSBU3bj0RecOV4Q/3uJl7bQxwbMDzfdgN319EZonIRhH5UUT+FlB2DxFZ6Jb7m4jcFu7+DZ9QVfvLoD9gKXCSe9wcp1U9yj3/ALjaPb4SWAy0BmoBrwHPutdaAgpUCsj3RWAITkOhGtAlTPntgK3AyUBlHHPOYqBKsAxh0n8G3AdUBY4DtgDPBcsF1HCvtQ1IOwfo5R4/AMwA6gO1gdeBce61E4B84G63nOpAd2ANcKCb97NuWfvGkN9I9557ANuAeu71Se59NwVygWPccpsCG9z4OW6dbQAahaiXHOArYBhQxX1uS4BTA+rtUve4FnBUhGd5OfBxwLm691bHvf+dwHtuGXWBhcBlbtwGwHluHdUGXgb+E5BXqecL1ARWAFe4z+1QYD1woHt9NXCse1wPODTdv6FM/0u7APbn8wN1lP5fwCZgGfAIUN29VvyDdH/UfQLS7Qfsdn+YoRTFM8AUoFmU8u8Epgec5wC/AScEyxAibQtXedYMCHuBEErfPX8OGOYet8X5CNQABOfD0yYgn6OBX93jE4BdQLWA61Nxlbh7vq9b1r4e89seVF+/A0e5978dOCTE/Q7A/dAGhL1dpGCDwo8ElgeFDQKedI8/AkYADYPihHqWl1NW6XcOOP8KGBBwfi/wQJhn1hH4I+C81PMFLgT+F5TmMWC4e7wcuBaok+7fTrb8mXknMzlbVfdQ1X1UtY+qbg8RpwnOR6GIZTgKf68wed6Bo/y+dE03V4aJVypfVS3Eaek19SB3ExwFsjVIrnC8AFzkHl+M0+LcBjTCUf5fuWaTTcBbbngR61R1R1DZKwLOA4+95LdBVfMDzrfhtLgb4vSMfgkh/z7ABUV5uvl2AfYOE7dJUNzBlDyvq3B6WYtEZI6InBEij0isDTjeHuK8FoCI1BCRx1xT1Wacj80eEn5G1T7AkUFy/x1o7F4/D6ens8w16x0do9xGjNgAVvayCucHWURRK3stIRS0qq4BrgEQkS7AuyLykaouDpHvwUUnIiI4ZqbfPMi0GqgnIjUDFH8LnJZoKN4BGopIRxzl398NX4+jqA5UZ1wjFMF5rgaaBZw3Dzj2kl841gM7gDbAt0HXVuC09K/xkM8KnJ5F21AXVfVn4CIRyQHOBV4RkQaEr7t4uRWnV3ikqq5x6/4bnAYBIcpbAXyoqieHkXsOcJaIVAb6AtMpXfeGz1hLP3t5EegvzsBpLWAs8JLbWl0HFOLYdAEQkQuKBuyAP3B+3KGmTk4HTheRbu4P+VYcG/Gn0QRS1WXAXGCEiFRxPy5nRoifD7wCTMCxtc9ywwuBx4H7RWRPV/6mInJqhOKnA1eIyAEiUgPHdl5UTjz5BaadCtznDsTmuoOpVXHMU2eKyKlueDV3ULhZiKy+BDaLM/hc3Y1/kIgc7spziYg0csvb5KYpIMSzTJDaOB/ATSJSHxgedH1tUFlvAO1E5FIRqez+He7WcxVx1gzUVdXdwGZCv1OGj5jSz16m4gxWfgT8itMavRHANZGMAT5xu+RH4Uz9/EJE/sIZ9LtJVX8NzlRVfwQuAR7CaeWeiTOFdJdHuS7GsV9vxFEoz0SJ/wJwEvBykHllAM4A8ueuGeJdnBZqSFT1TWAiMNtN95l7aWc8+QVxGzAfZ6B5I84Aco6qrgDOwjHTrMNpFd9OiN+lOmsTzsSxof+KU7dP4Ay0gjMQvcB9Pg/iDGjvCPMsE+EBnIHv9cDnOGauQB4Ezndn9kxU1S3AKUAvnF7gGkoG0AEuBZa6dXodzrtjJBFRtU1UDCMYcaaZfg9UDfqYGEaFxlr6huEiIue4Jod6OK3R103hG5mGKX3DKOFaHDPLLzi25evTK45h+I+ZdwzDMLIIa+kbhmFkEeVunn7Dhg21ZcuW6RbDMAyjQvHVV1+tV9VG0eKVO6XfsmVL5s6dGz2iYRiGUYyIRFq9Xown846IdHe94y0WkYEhrl/uet2b5/5dHXCtICB8hvdbMAzDMPwmakvf9akxCccD4EpgjojMUNWFQVFfUtW+IbLYrqodExfVMAzDSBQvLf0jgMWqusRdVTkNZxWhYRiGUcHwovSbUtrj4EpCe0w8T0S+E5FXRCTQYVI1cTZ3+FxEzg5VgIj0duPMXbdunXfpDcMwjJjwovQlRFjw5P7XgZaq2gHHJ8nTAddaqGoejk+VB0SkTZnMVKeoap6q5jVqFHXw2TAMw4gTL0p/JaVdnTbDcZxUjKpuUNUix1SPA4cFXFvl/l+Cs8FCpwTkNQzDMBLAi9KfA7R1XfBWwfGWV2oWjogEbvrQE/jBDa/nupBFRBoCnXG2XjMMwzDSQFSl7zqc6ouzjdsPOFvhLRCRkSLS043WT5zdlL4F+uFsxwZwADDXDZ8NjA8x68cwDCNl/PzT9/R/6m7+2rIl3aKkhXLneycvL09tcZZheGf4xOE03K3ceOvIdItSIbjylYnMbHAcly55kwlXDUq3OL4hIl+546cRMd87huGRcXf1p/Hsedz89D3pFgWAJT/9xMX/eZTHDj6Hme33Tbc4FYYtVaoBsKNKuXNIkBJM6RsZzZBHR3LYu2/w7OR7E85rd+2aAHzWpF3CeYVi7OhbGT4xePfB8Dw781ner+vsI76qyp5JkcnIPLLzU2dkDa/u14VNUp8Vv7+XblGiMuWY89ghNRjhMX5BrrXZ7r/vLqpXqsp1/eIw00io2eiZjyl9wygn7JAa6RahwnF3J2e953UxpSpf45ipxpoKWchz/3qIh+4elm4xDGDCqFu55am70y2Gr9z14DCum3Z/usVIOm3e/x+XvTYp3WLEjCn9MAx4fAwDHh+TbjGSwqMtmjDmiHPTLUbMjB51G/2f8j6IOmF4fzZJ/SRKlDhvH34YL+xzKmPuHhBzWomzwTp+7O2MvO/O+BJ7YHKHc/nPXicmLf/ywlapzdv1OqdbjJgx804Ynt73dMDZHTvT+KVSGU8YFYLnOp/OJqnHQZPu4aob7oga/8VjT0qBVInxZ6VaABTm5qaszIlH9aJQcom1r3fkrP/QcttaXjrr2qTIlWqy1chjLX2jwrCF2iHDr57+IIMeG10mfFVuKL+AqWXw5FE8OH5ousUoRaHE94FZVqklH9Y50mdpUk92Dt+WYErfKLeMGX0rtz8xNmq8Nxodz5PtzkiBRGUZed9Qxo+4OeS10eMHMXW/M3nv4IrZs8pU0tnC3/DXRr5c9m0aJTClX65p+f6n3PB84vPLE2X4g8NoPHseY8ZGN6n4yTPHnMmzbXrElOb2f42j8ex5jB15S5KkKs0jnc7ntc6h7dcFVZ0W9foqdePKe+S9Qxk3+ra4ZasIdHnnFS7+z6PpFiMmvps3h2defDyutOd8/j49l6TXsGRKvxyzQ2rwapNujB4aakOy5PPx27MYO7Qv81u3AGB9k4YpLf9P2SPmNF83awnA1kaxp42X5bn7JCXfRw49nwc7X5KUvEMxcMoYbnomttXGT0+5n3HjYx+ELmJx5X2LF5glwvJlSzzHDTTv9H/mHu6fHL03Gcg1637jjsaHx5SmiJ9y079y2pR+BeDhbldHjxSGfs9OYOS98dmUn960gIndrmZL5epxl19RGfbwXdw/rsx20BnNU21P56Xmp8SUZkqrFjx45EWe48/874tc9J9Hef7pR2IVLyLD5/xfXOlebH4Kd+8XW29yWaWWcZVVXjCln+FMb3Yyjxx6flxp5zVoDcCflUIPoMbKuLtuZtxdoe3f5YkxY29nyoFn836H/cLESO1Q4NjRtzHigRBTLGNwljji/jsZP/Z2H6VyiHUm2OvbVjO77tG8XdtfE8cvdfbyNb9MxqZsVlCefOx+FlbawRE19uKCi65MtzieePD4ywEo734NC6tUBmBt1QZplsRh8jEXsEuq4sUrz2cfvc+U3xdw8M+ruGXQuOLwRzueR67mE0/f5YN3Z7Ltzz+g/oHFYXc8PpZDcmtCq+NjyqvQdX1QUB5cIJQHGdKAtfQrKF/UUZ5tfRo3Nj6Uy195KN3iZBmpHYjb5exD5IkZCz/hzQbH8kGHsk7hCiS+Nl5fqcyVAQp/1L1DeGbfHkzfy58eYCgmTRzDA/ffFUMK7wrc66K2J5+eyKjHHHv/yCfG8dGsf0eM//vatfR+dRLzvpvnWZZ0YC39FDNm2E3UrlmXfgMS832+pVq14uO3GhybqFgRydZFLMF8fkDbdIsQlULXCdvunFxOffM5qhbuYsbpifUE1+eU3rd6t+uS+I/KtRLKNxKjDnYWR4YyBo6YOJx+l9xEvfqxrba+/sX7WFurDpU9xh/U4jgA/vb9VzzS5jTeKFjGlxHid1i4Gup35tdVC5jVoWNMsqUSTy19EekuIj+KyGIRKdNDFJHLRWSdiMxz/64OuHaZiPzs/l3mp/AVkYdOvIJ3D2qVbjHCMnbELQx5dCQzXn2WbF/GEvyx+7r6ITHmULr+7h5xC6P+OThCjERyL8u31Q7iyxqHhr2+bIn3GS/xMvnB0Xw8+y3f8hs4ZTSPHnwOw2ZOjTntvxt35dNaUfcYKcP2rZsBWBf08QvHjpwqMZeRSqK29EUkF5gEnIyzSfocEZkRYtvDl1S1b1Da+sBwIA/nN/SVm/YPX6RPMVe+PJFPGnTgx64nJJTPlzX83xv+lWlTOb9XfC260fcMZFPDuvzzykG83OUkVuc04TndwU6pFj1xRSGN9tuikp/vciq/5+xF8rzelKAePifPvPYEObvz2dKgDuN7l8zwmvn6NHqc2QuAL/73HhD/2MZdHc7goJ0LeTfuHErzZ23HE+mGWjV9yjH78NLSPwJYrKpLVHUXMA04y2P+pwKzVHWjq+hnAd3jEzX9zGx4XFxzx/1m7Kj+bK1c2s77Qs3CuPN7+PBePNfqNADWiDMLIrTCL2vo6THzaRrPLt82zGj0eDP0PazIbREy/rJKLek86zU+/987MZXze07yZ5jE4oRNJYdJR17IU21Lr2a+stb+xbOFZv/v7YRl+r5q+7jTHvLeW9z+r3FR48VjgnylSbc4UnmhfPeQvSj9psCKgPOVblgw54nIdyLyiog0jyWtiPQWkbkiMnfdunUeRc9eJna5jM9rHlYqLJ5u6xgfFn0VmTxG/XMwjWfPY/Q9pa1/Nzx/L4MmjyoVNnzicP7+b3/naRexOTd2O/PX1cKbbSaEWdn7S6XWfPKBX+3XeEl8tCWcH55HDzkv4bzjaQy89+7/lUq3Nqcx01qVVc7Lasa/U1hhnL2+7VKDkU/+M+5yywtelH6oGgp+214HWqpqB+Bd4OkY0qKqU1Q1T1XzGjXyZjdLJ/2emcBtU6O3PlLNnQ/fFVP8hzwv+or+I/m5tfMtf7PT4YweVeI64NUm3XhyvzNLxX3s4HN4b49jyuRx+LszOP/1ssvbh08c7nkbwRW5Tnvjj7qlu/8Tu/zDU/pksDa3EcfMei3p5fR54V6eant60sspwf8W7bxvIw2VlhC8PmC3eB2ehZXV49cxj7Qs/55bo+FF6a8EmgecNwNWBUZQ1Q2qutM9fRw4zGtaPxl5352MGBjbHjrxML35ycXmkCJG3z0wJt80gx8dFfH6e/99I2a5Hj/w7JjTjLp3SPHxkEcTm1EEzo/x4S7xuQ5YkduCj2sdzvkzHmfsiJIW9mMHn8NjB58TU17xtObGj7iJa156wPfpSn9JbZZUal18/uP8+f4W4PLa3kkwV4Soxnj9+AP8Xsu/aZ7Tn5tSfLwziYOn69atiSm+1+p5fdlnNJ49j+cXz45dqATwovTnAG1FpJWIVAF6ATMCI4jI3gGnPYEf3OO3gVNEpJ6I1ANOccN8Z8TA63ik03l8eGJsi0X84uEjevHQ0ReHvT52VH/GDr2x+PzNdoeFjRuK0YP6xBT/xuf+WWamSCgmHXpB8fGsth19a7tF+6hF4uPah/PJ4WVNLrHMM4/nPh447gpe3/MENjYqcZC2vHWTsPFzcuKrrQk/vF98XOD+BDUnsSUzEqRqNlauU3w8boy/jvI0gZfks5qxmyEDKQiop2FNSk+hnfjQGAYmYeOjR3Vr3Gl/WvYz7d7/gHe++rDMtceWrwXg1hX14s4/HqK+aaqaD/TFUdY/ANNVdYGIjBSRnm60fiKyQES+BfoBl7tpNwKjcD4cc4CRbpjv5FZxBh6XV24GwOjxd3DUrP8wZthNySguZiZ2uYyJ3a6KO/0Hx3eJKf7LTU9i0mF/iynNDqkat6/1YLbWjL6gaMT98fkEeqdZ3ZD24utevC+u/ILRgF7Cy00T6M6HUY7z6pe0+n/Ldd7XRzrF7iojsMUdPFsn0AnczlphZmF5UN4SUcN7b/KPmTCEd2b+x3P8cLzRqGRNyiYprSz/eWA3ntr3dO59MP4GRyh2VPJuOgpm2g+fsFn24LHffy1zbZumbuOcQDw1L1R1pqq2U9U2qjrGDRumqjPc40GqeqCqHqKqJ6rqooC0U1V1X/fvyeTcRlk+79CepZVasrHV3lHjjhk/gFETQreK73pwGGOH3Mjo8f62lv6UOjG1hn+oUnaFpd/8KfG5AIb4LCKPdjw/Lmdw4Tby+E/jrjHnNfLess89v1LFWKieSIvbK9OO9Kfn/FDeBfyjesvi8/f2OIaT33reU9rdUqX4Ix+uUbI6pwm73BlnEzqUjCH1ff5eT66bR04ew8lvPU/vlx7wJFMgTz7/BFs2bSo+V4+mxYUFzaNHSgIZvyL3+ZbdieaR/iHXS2DuOEexDx5U4l52codzaVj4O+tzYp8tMHDKaPZYsoaB4x8uc22b1GLq/mdSa2R/Cnbu5M4xpWezrN+yBuo0Kz6Pdwm9H/hh4h41YTDkle15bA/TI/i62kE+lBqdRw4tK9MrTby17nNzM/7nw+qc8OatRJlf9cBS55MmjmFJ4zphYsdH0bTM5559FJo5LpxDrWF4ZD9nAHz+njAl+GKUH8CgJnnMnv0CuJMTCij7YfqkWkd6z3wCqidm3vKDjHlrc6Qg5jSnvvkcHX5bDu5GHROPcmzywW2/eBT+8IeG89RB59Ct0WcR40089jKaFPxWvGBn7Lg7IH83bVodGDFdrPzjtUlcXW1fjutxatg4wXbhIlbmJt4iya8auossYVpF6pOZyQ/uHXUHa5rUg9alB+9zErTDJ8ququV75WesPNW+Q/Hsq1CMuW8oeDSDjX3wLqpKDhzsWKD/tWeJKSgZbkVW1ChxCbE0jOvlN6t1SELJsZMxSr+IcIorFN9WO4hv2ySnRbm9mvOD/KNajahxA/dyLfrwPPjb3LjKHT1+EB92PAiCWlHv1OvM/u/+K6LS97KK89pp99Ni6RqGDMzELeND89Gh+/NFBHcGqWTw5FEsaN6UVuvXUyXX20dnY72SGTODHx0F+zvmjxkdSsxkf5vxGNSOdf9bf+1L63Mir/zdVt2747mJHUrPZNuR4z2tnxQGuL/eLeXjI51xSj8UIwf3YdhYx3wyasJgdtSqxpjrhyWlrM6zXuXoFT+lzX3pm4flhfVxnqgNeEHjpnxX9SBq77mZrY+NRgoLYf+eYeN/8f4HULSCOYzv9x1VK/H04w9w2TWx+9kfN6wff+3dMKIMfvBn5dBL/nOCBvgufW0StXfs4JGLbwXgm31bxl3myHuHQoh9EKa6ax6+aAHH/FXSMFib0zhsXtObn1ySfv8Se3dgD+6jGBR+pCmb33/3DeVxReruGMyjjWfP446f3oZ24RtIwfxQZf9S5/3+PYnpe3SGav67XEmUjFH6EsEc8MjJvSlS8ZNcu7L/E7scfqnUhl9ateEfi2PbyWfw5FHFP2iAxYvmQdN47H/J+8F9V9XpFe2kClODNiL/pm7p6XN/1qjOK79+VsYkEswL+3Rn+ZY5xOOJ78ETY/M1dNeDw1jepCFTL+gXU7qtOdF7awCz6nUGoGh0JnjVtFd6T7ufGR42vtlYNXmujb2wIbfEZNJ72v20WrqG+fu3hDi3P0ymN9f8AKXvZeeruXuVjC3E7mgPR+GXUyrGNIWYSPzVGXOPH9t8xKZ8P29Veu/MUKtlxw65Mal+bmIxjQUTPO6xtmZdNtYpUZYFlcJ/lOfW8t/E9umeZf29TO5wLjMbHhdzXpHszEV4XTHshRl7hd5o/dYnx5c6X1Q53M5eyaWox7ghp2TP5Bl7ncjUI3owv3ZsO2nFQkECs6q8mC6zhQxU+onz0OEX+peZj82XXXtGd/a2M4bl6MHskOTthftzs/Dmh2SUG8nc4RdVA8w7sawYjlf9vNMiteMK577+REzxt0hd1sUx6aGIaPUyvU3qFl5uT2Bufnkna5T+6LuzY5PrSDNtdu1Vn3+//HTY636ytnLpQblwc+srMpWqhHc9fdqbz4S9tiK3uWM3fmJsTOXtiGEHLT+Ix4lfIkRrH22T5G3aEkyiK4fj4doFS7l2wdKkl5MxSl8l8q08fESvFEmSPLbVSMy//WMHncObBUlZEF2GwBlJmcqA/UKbYQC+8TA97xl3qrBXtiSweK6844ffp4rMj1t38N/fN/Hf3zdFj5wgGTOQW4QAo0fdxrqjkruFYDp4Zt/EPSi+vmfsNu1gYtiNNOGyksGEkXfw035NaP3DUnA3a/eTB8YOgqN9NBFmAf/avyfVdFvS8i9vNv3gsbnjv1wUJqb/ZJzSBzx5ebxz0l3QPnavlIHcPeJ2CnU3g+4qu3R7d+XYFhctrHJAQrJ4pTwtegpk7PB+cEJie7l65d5jnbUQ7ff4IUrM+HirU/ybhlQENtRPzqyhHeJtlpSRGBmp9L3weIIKH+D+4/4OQP7g6+Hka0tde3EfZ45v+WpflF8mpkjhB7IrgUHvSCyt2ix6pArMqx7dVJQn1knD6JGyhIyx6ecQuxsG38pO83L88srqmplrg45EsPdHI/2U1x5uOshAbZXMJR6hEc3AaoyA1w3T41nUkkp25WTutDzDCEfGaavNETYuv/yVh5JTaBQbzvkzyvjtM8oBgT7nDSNbyDilH4m3GqR+Ro+gfFz7iJSXaxiGEYqMUfp+7fgUFzZaaxhGBcGT0heR7iLyo4gsFpGwS1tF5HwRURHJc89bish2EZnn/k32S/AyaGHSso5GQZ3wU9jmVu+YQkkMwzAiE3XKpjjuKycBJwMrgTkiMkNVFwbFq42zP+4XQVn8oqoZrfmeO+yUdItgGIbhCS8t/SOAxaq6RFV3AdOAs0LEGwXcA+zwUb4KwV+kzieIYRhGInhR+k2BFQHnK92wYkSkE9BcVd8Ikb6ViHwjIh+KSMiRVBHpLSJzRWTuunXrvMoenEl86XwgnfvXGoZhxIIXpR9KmxZPhheRHOB+4NYQ8VYDLVS1E3AL8IKIlNn5WFWnqGqequY1atTIm+SGYRhGzHhR+iuBQH+9zYBVAee1gYOAD0RkKXAUMENE8lR1p6puAFDVr4BfgHZ+CG4YhmHEjhelPwdoKyKtRKQK0AuYUXRRVf9U1Yaq2lJVWwKfAz1Vda6INHIHghGR1kBbYInvd2EYhmF4IqoxWlXzRaQv8DaQC0xV1QUiMhKYq6ozIiQ/DhgpIvlAAXCdqibHobvNlTcMw4iKpxFIVZ0JzAwKGxYm7gkBx68CryYgn2EYhuEjGbMi1zAMw4hOxih9UbPvGIZhRCNjlL5hGIYRHVP6hmEYWYQpfcMwjCwiY5R+6vfLMgzDqHhkjNI3DMMwomNK3zAMI4vIGKWvOTZl0zAMIxoZo/RzMudWDMMwkkbGaMoCtaFcwzCMaGSM0jcMwzCiY0rfMAwji8gYpW/juIZhGNHJGKWvORlzK4ZhGEnDNKVhGEYW4Unpi0h3EflRRBaLyMAI8c4XERWRvICwQW66H0XkVD+EDo3N3jEMw4hG1J2z3D1uJwEn42ySPkdEZqjqwqB4tYF+wBcBYe1x9tQ9EGgCvCsi7VS1wL9bcMvausPvLA3DMDIOLy39I4DFqrpEVXcB04CzQsQbBdwDBGrfs4BpqrpTVX8FFrv5+Y7mmqXKMAwjGl40ZVNgRcD5SjesGBHpBDRX1TdiTeum7y0ic0Vk7rp16zwJXobC/PjSGYZhZBFelH6oyZDFBnQRyQHuB26NNW1xgOoUVc1T1bxGjRp5EMkwDMOIh6g2fZzWefOA82bAqoDz2sBBwAciAtAYmCEiPT2k9Q1RM+8YhmFEw4umnAO0FZFWIlIFZ2B2RtFFVf1TVRuqaktVbQl8DvRU1bluvF4iUlVEWgFtgS99vwvDMAzDE1Fb+qqaLyIPo62vAAAgAElEQVR9gbeBXGCqqi4QkZHAXFWdESHtAhGZDiwE8oEbkjFzxzAMw/CGF/MOqjoTmBkUNixM3BOCzscAY+KUzzvmhsEwDCMqZgg3DMPIIkzpG4ZhZBGZo/TNvGMYhhGVjFH6hZKbbhEMwzDKPRmj9A3DMIzoZIzSzxGbCWoYhhGNjFH6hmEYRnQyRumL2fQNwzCikjFK3zAMw4iOKX3DMIwsImOUvkrG3IphGEbSME1pGIaRRZjSNwzDyCIyRunnYPP0DcMwopExSt8wDMOIjil9wzCMLCJjlL45XDMMw4iOJ6UvIt1F5EcRWSwiA0Ncv05E5ovIPBH5WETau+EtRWS7Gz5PRCb7fQOGYRiGd6JulyiOf4NJwMnASmCOiMxQ1YUB0V5Q1clu/J7AfUB399ovqtrRX7FDoIVJL8IwDKOi46WlfwSwWFWXqOouYBpwVmAEVd0ccFoTUP9ENAzDMPzCi9JvCqwIOF/phpVCRG4QkV+Ae4B+AZdaicg3IvKhiBwbqgAR6S0ic0Vk7rp162IQ3zAMw4gFL0o/1EaEZVryqjpJVdsAA4ChbvBqoIWqdgJuAV4QkToh0k5R1TxVzWvUqJF36QPz2L07rnSGYRjZhBelvxJoHnDeDFgVIf404GwAVd2pqhvc46+AX4B28YkamSGjJyUjW8MwjIzCi9KfA7QVkVYiUgXoBcwIjCAibQNOTwd+dsMbuQPBiEhroC2wxA/BDcMwjNiJOntHVfNFpC/wNpALTFXVBSIyEpirqjOAviJyErAb+AO4zE1+HDBSRPKBAuA6Vd2YjBsxDMMwohNV6QOo6kxgZlDYsIDjm8KkexV4NREBDcMwDP/ImBW5hmEYRnRM6RuGYWQRpvQNwzCyCFP6hmEYWYQpfcMwjCzClL5hGEYWYUrfMAwjizClbxiGkUWY0jcMw8giTOkbhmFkEab0DcMwsghT+oZhGFmEKX3DMIwswpS+YRhGFmFK3zAMI4swpW8YhpFFeFL6ItJdRH4UkcUiMjDE9etEZL6IzBORj0WkfcC1QW66H0XkVD+FNwzDMGIjqtJ397idBJwGtAcuClTqLi+o6sGq2hG4B7jPTdseZ0/dA4HuwCNFe+YahmEYqcdLS/8IYLGqLlHVXcA04KzACKq6OeC0JqDu8VnANFXdqaq/Aovd/AzDMIw04GWP3KbAioDzlcCRwZFE5AbgFqAK0DUg7edBaZuGSNsb6A3QokULL3IbhmEYceClpS8hwrRMgOokVW0DDACGxph2iqrmqWpeo0aNPIhkGIZhxIMXpb8SaB5w3gxYFSH+NODsONMahmEYScSL0p8DtBWRViJSBWdgdkZgBBFpG3B6OvCzezwD6CUiVUWkFdAW+DJxsQ3DMIx4iGrTV9V8EekLvA3kAlNVdYGIjATmquoMoK+InATsBv4ALnPTLhCR6cBCIB+4QVULknQvhmEYRhS8DOSiqjOBmUFhwwKOb4qQdgwwJl4BDcMwDP+wFbmGYVQI6hduSLcIGYEpfcMwKgStd66IHsmIiil9wzCMLMKUvmEkEdHCdItgGKUwpW+klRP+/Dx6pApMs8KV9Fr+TrrFyFqaFPyWbhEiUp8/Ul6mKX2jFI0Kf09pefW3bU1peenggcvuSLcIGYFomcX8Ufn6pNOTIIl/LDzxxJSXaUrfKEW3lfNSWl7zJWtSWl6qCeWHxPCPg3cuSEu5DT00jvbb/VMKJIkdU/pGKTTFWkoL/VurV9678kZpum/4X8J5dFizsvh4zYkdS10bsuD/Es4/EVrv+rNM2IBqzoegvm5g/lHpcS5pSt8ohSjc9MnztMxfmm5RYqZhfurto5lGFd2ZusIk8RZG1Z27w167se+QhPMPT3yy719vbwBq6jYaVa/vp0CeMaVvlGHQ0Amc/t2cdIthADU0tWMeZ6z5OKXlJYuK2OurUymHxlUqJ70cU/pGKYrMO3feOoZDdnyfkjJ7/v6B57g5ruum2lq262xUMOIYmPXKzSt+TVreyWJRl4P5+pjgTQn9x5Q+cNTWr7j5oyfTLYZnrvzx9ZSUkyr7fvWduzzFu2DlLC5c/h4AB2/9OUpsb7TOX+JLPuGoVZBYS737Gv+ntHb98zPq6Cbf840ZH8w74fjHFX2Lj2voXwnnN2rp+wnnAXBks4PYQzfSt17Ze88RISeJdVJcTtJLqADUyN/FwOEPxp2+nm4sPvZjcCoaY6+709f82uT/UnwsyWt8hUSA3MIUFxpA0+3rk5p/t6++jhqnqu4IGb7mxI48cvGtfovEwQuW0GnLj77nG4rr573iW17xqsPaPij9a664pdS5lzc2J0Ss+rXqsahrVy479IyEZYoXU/o+sPfutcXHlQtKz0bp/97TqRYnZi5dlBoFEI4DNsen9C/59U3ytkeeYtqsIL3+WoYMnlDq/OitX6VJkvQwvP/odIuQNtLXlImMKX0Sb91WLQw/g2DA6PvTNpfYK9fdODit5V/Tb2Bc6f555SDarY1/nv91370ad1qAPQvXRo8UxEHL4v8IHbEteq8hXRy2/dsyYYPn/Tdl5XfYGd/4U59Fr3PHopnRI4agbnkwkcVBVij9Gz97IeL1RG3XjbdGHlSstyvx7qXfJKrw0kIUe+exC36ImkWgKeWum0YlJE6rHbHPEOlUr3n0SGGotTv50ykbr4lv2muodlO//sNjThNrKYfU2JOauoVjF0Uf4wnlB+mY/Tpyy/XxNXr6/7aoTFiXv+bGlVcqyQqlP2TwPekWodyRqMIrNwR8CO68NfRePQfuKvkYSBnFkdpO+LkXXxV32sabNvsiw551G4YMn7xhHsNuCW+OCazHZCBxbKp34aW9+aXrsdx58wgPsaM/a6+raPfKDz0WNKXLuWXC4ukRJhNPSl9EuovIjyKyWETK9MVF5BYRWSgi34nIeyKyT8C1AhGZ5/7NCE5bEbj2u9ciXt/3p+XFxxqiNdpuRXpcDZy9dnbSy7hoWWLOxNS1rQUOJoePHJ+CrpFf0ro/b9mHceVRHrjhqLPo+8VLCeWxV+Earuo7gGAFOGDOK5x9/uUAvJoTfbOSmz57MSE5QnFs0lvJ0bv0lQvzOfP3D7hw5azisL0KS/9++/36Lvc1bBU1r6NyatC4cDUX7q5gSl9EcoFJwGlAe+AiEQmeTPoNkKeqHYBXgMCm9XZV7ej+9fRJ7qhcttj7EuxoNv0qW0PPrigieLAumNE3DOfc1e95lscvav+1PWR43y8cc9cZ6z7khq+mJ1TGXr+uSih9ETlhXBA3KFzn7YPgcvTWEsVxwubPOefLj0pdv/eK+MYPvFBTt8SVzuunrM3++zN04Dja+DzNtOfa2fS/o6SF3/n4biHjBTo8GzT4bipp+LGseBizX2df8yuibO+uhFwpuxjq8QtvpvnW/OLzb7t1L3V98JW30fGwo0JnWFjyHtetVI153U6jU4NmMUqcXLy09I8AFqvqElXdBUwDzgqMoKqzVXWbe/o5kPa7rLMuPYt3gmfvlEeGDnS+yU/87SbuvG1sqWupnrIZjb9/+iZNtnvfJq9qfsmPddpZ1zFk8ISUOT2rFsGFwdUL/1v8sU2UUSF07QG7ytqXyzvBr1rb/Q6KGN+v59iiYHn0SHFSEMJEVZjERWjx4EXpNwUCpxysdMPCcRXwZsB5NRGZKyKfi8jZoRKISG83ztx169Z5EMlfog3kFm7xbkutnB9a6Ycy+wRz/bevcu38f3suKxg//Kak2uFaEZXC2HMH33lflJSJ/aBSdbujbxhe/LFNlK49ytqN08Hh274rOUmwImurP+MVkXHelXMXptaTbHnDi9IP9ThD/tJE5BIgDwi0d7RQ1TzgYuABEWlTJjPVKaqap6p5jRo18iBSdKQwtLng8G3fxJzXneMmJypOWKoVd5Bg+M2jGNHPy4BUaJZ3PbLUed2lsc8uCaf0b//fc/GI5Jl9/kpskVSDwtQ3FoLpuPknnt2+NGUrptPNkAZt+VuA7fv6efHPCLthnjfzZ59vX+XGbxJb8FU5J4dqGtr0GckUFI02+Uso9NFrbLLwovRXAoHzzJoBZQy5InISMAToqVrS5FTVVe7/JcAHQKcE5E2I+oUbeP30Kzh7Tekl1X6aNMpTR65Bg8Yxp6m6Mz9k+K3D/hk6QaKraYuSx90Fdr5S++6INP89OU+l1frSH6rLqzbh5B5np+0lOHbLl3Gl8zqFr7KWfjfyjj6OajtKXGgM718yI+yahf+JnmFA7/fmW7w1dobdPIohEWYYJYom1GUJ/eBrVq4GwIG7y4cJzsvzngO0FZFWIlIF6AWUmoUjIp2Ax3AU/u8B4fVEpKp73BDoDCz0S/hIBOuiPu89Qa9ZzsyHyRfdEiJFcoln1589dCNXL0zdAheA3DA9pFCc9MenHH/8WWGvX7jiHeoXRrHHe6iXyDG812u3TZ95juuFaw4+udT5yT1CWi9TRrzfmmN3ePPseMp30acz9vn6ZW7+9AVG3XBXmWt1yomTPPX0zkWPkxPiA1EYoHja7OX4yz+xU2f67JzH4+2PLBMf4IxGdem/z15Ry/OLStEiqGq+iPQF3gZygamqukBERgJzVXUGjjmnFvCyOF/v5e5MnQOAx0SkEOcDM15VU6L0gxk2+uF0FJswlfJDt7yjcfr6xKcmViso7QjtxjkvsXbv+kxv5ii7587t41yYHdpG+uA/7qD9ex5nLXkY84gUw8tHtfEmf5XOgYd2CnvvFYXzVr3L3y+/zVPc/gNGc/fseRy/+QugY8g4w8KslQDo9/2HbN65nYl5F8YjatyEMtkcs3ke79c9mmZNkzfnJO+QI4qPh3W/PGy8Jw6KPv3TTzz17FR1pqq2U9U2qjrGDRvmKnxU9SRV3St4aqaqfqqqB6vqIe7/fyXvVuLngJ/Dj+bvk7/MlzLCDeQes7ns8nU/+NcFN7FrZ+wDu4E2/SM/Kz3+MeSOcUy89PawaQNnp1zx0xsxl50s9vrLGSSsFWIKa3mbrZQIqRiU/nrfBtx/aI+40vbtN4zBt4+Lu+wmm/1zezDpuAt4p2E+bds5M4bCOb2L1QtrgdtTrhImv/JAxq7IzfXoorRV/q8MvcN5Ecf/8gkn/Fnale2Zs9/ylM+Nnz3PlT++HtOUzW6bPuPMLcnfNCGQWGb4DBn9UEx5B85OGXft0JjSFpl5Gheu5vjN8bsTPuinpWXCzq7SmGu/e40RN8Y/SB4LiQwGlneaNG9Ok+bxu5IIpJ2HNR4dd8wH4JqF/2XsKZf5Ui5AvfoN6XBwXvH5pYtjd51cGOE5l+e9kTNW6cfD5VffUMYd6tCxj3pKO2TwBMZedyctV2/hzAibgnTa8R091jsLhnILC7jo0qvjljca+bvK+qk/9s/UenkM9bNoWrAyRKjDYRt/4qWzrosh/9I/r6F3jC8T5/RzejHippGe84yFi5d5axREot2uyK3JqxbNiNuhWLz889dPuDmKz6pEGRrCFHTT3Je5LmAFfJ+tcP13rzHqhuHUa9DAp5I92PTj9GuvMYyJpYuoNv2KSiH+Vv618//NYwefEzVev0HOi9w4jK233s6tngdL221f6lm+Y7c4dtbuGz+mwea/gI4U7vZ3xWRSEH/aHfG2rPbcsgXqxF/ufZcP5IUE7fo95y/ihIL5PHxkaFv3mOuHlQm76ZPneLDzJcXn+2xZx9YqVROSI5BLrrwBgAdSPGYx6Hbn91N1whB2VKtEzxtH4NcyfsHfiVWhBnJLyiq/vb2sb+l7naI1ot8IWhQkbt+P5VXo/Ol3nm3OVQqcD8lT5/WN4mqgPHc8U8817Y4PGX7QztTNN+h3+0iGDozP1l2/cAM3fvki1zSMdZu91L4HsbohHnT7mIjmuBPCmADb7l5cJuzShV/R5a857KXRfeAUKet4ZtsBNGzYiBwtoNem2NcDpYqMVfqF27ZFj1TOGTDmgZSVVWQ79VsZ5PrS40qegupwxOEhw49bWFbpx+MFMtkoMGTA3XQ+IbS/nGikQvXf9sV0bl80x9c8p511HWtOLDuD6P46jbn3ty9Khd1y4528cuY1Yf07hSIW884++Us5fsuXDK9Uh+rVarCq62GML5rZVg7JSPPOxcve4s6xjzIpga5p9w3/o923P8GJsQ1mRiPb2tmV8dHElKIec13dRJXNZfdAuGHOy6zbuz4vNT8lofyP3PY1HZcs47GDopsLY8MPD/WhOXLb14SbpunkE/7Nvm3g2LDXrlo0g23Vq0bMOxbyjupC3lFdQl7ruehrXt9fOaXLaZ7z8/J7zUF5qWdvz3mmm4xU+n5QffduBo8KVvglr8B+u3/ix8rtEiskyZsgF0h8c/yj0efrV9AcgRAtrXDkaAGFkgtA59U/ML1Z5PnRJ//xCbPqxed18bLF/0dhTk5M8hVRTzfyQ9euPLe07MDp0AHOIPFLURoT0VqJ/z39SgAe88NeHvc75M8XNNE3ONR4RbIY3ncYkbd1iY0U7GGeFEzpx8mRyxdzhPxCrc3b4lIufhPKZ46GcJHgh0O1Ybd6Xwbfc8EXzGyfxy6pzFppTN8vpzF0wHimBym8dvN/YsUxDWk9/1e4AJ4994a45bv7miFxp810KvvsDjmbKWfOMz1jSj8mnKd8w9cvh92lKdXU0i38JbXTLUZYRvQbwQjgsHdDL9bKce3kg0c8xGCAUy9KmWxF9J3zEjtrVuXx9ul1o5AKzvviIzjpzHSLUe4J5aph78JVHLh5KX6ZotJFhg7kRm/O1nN9wlQrTP6+owkRpTVx0lp/Bsi6fLuQ5gUrqL3y9+iR4yAxR1beSoiXoXeMo/oGx0VDLIN9iXLC5s9j2uwnkFpVawBwxB/eZxnV1C0MGhJ5w59wxDubJZP4plsPnjun/A7QeiUzW/oeXtAzf/2CZ9r0oHp5V/pRiKRKi6qhsu7ijLUf8+/GXcPGHTpwHEMha1uBZ5z2dxat+IT9fl4O3U5KSZnTYliEFsyNt4+k9iN3k3f0sZ7it85fwomLv4eu3uJnO1JRDfYeyEyl74Usarkke6FIv/eecEo4MbxTu2AZ6ugmNssenvKvUugMSFcqLEDLdE79+XF2ODyPZw4vWZZfo3rNuPMqahV32vEdyTQFXN5nQOlyI8T99ORz4eT4N1+pWpCcSQHlh9LvpxdPnBX1u5Ch5h3vxKYQK+hTBpr/tJImBb/Rfr73/Wa9Mnj0wwyJ0YvpAdsd1wOFu6ObUw755kd6rP+IQwtrhY3jt+O0mrXjHyfp2uoQTvnjE477Ij0eOCO5AYmVLn855sP627b6lmfFIvMah9nb0k8jXj4dh+z4nsBWYlXdwU6pViZei19WQ2No/2Pk1cKDh9/nDJSedHpMsiaLE+cv5sj8hQz24NRtwIiSRWpfPeXPloPJpFv3s+gGEGPDuk3+Ehru2sQXNQ6Nq9yiD1/V3f61yqtkfAs/+8hIpZ+TNNNNbPl6XXrec+1sGv7xV/HUz34f/IvKVIHTLikus2nBKpZUal0m7aChExgEIaeN7n9IHu12/0yXXxZB19AbOKSLmyMs2EkF5XF17Seu+SWc36ZoHPn9EvI75HD45hxeibSLdQxkXjvXK9GbZq0bNAFg/83RvYWWJzJS6ZcHH+nXf/MKlbdsg64nRI07pVf/UueDR/izCvjcv10ea2OzwrHHH5thH2ixcaPnNH2+nk6lbTug62FJlCz13DLkbor2hBvgfjhECzjz949IdGwh3OydOn84K5dbb0z/HsWppttpZzLtrTfJO+XSdIsSE56Uvoh0Bx7E2TnrCVUdH3T9FuBqIB9YB1ypqsvca5cBRc7VR6vq0z7JHpZq20tm5FTWsu6FU8HwJO7jaZQwvP9oKo2+lR4nnu85zbBb09vLSCWrux4GJO/jNuSOcVQeN4BLLr0xaWWUZ07o7t2lQ3khqtIXkVxgEnAyzibpc0RkRtC2h98Aeaq6TUSuB+4BLhSR+sBwIA+np/iVm/YPv28kFH2+eZVKW/+CrkeUuVZ3o7Ob0v5r16RClJCUgw5JVM5e8z5v7nVMusUopvGvq+jQ+HvyvlkA7uzSIUPv9b2c2rVr4bN37ozljkF3p1sEIwa8tPSPABar6hIAEZkGnEXABueqOjsg/udAkaPvU4FZqrrRTTsL6A68mLjo4am00lHkw24ZFTbOkAF3UzC4D8PGPuI5384rF7F8n33I/WNzYgLGMeaQ/MVNofFrE/lTf/6GafvVo/K2slsWxsKAEQ8wAKD7JdGiJo2rFs3gr5rVqOgrM9PF1T/8l9z8wnLhviQb8aL0mwIrAs5XApFGBa8C3oyQtswQk4j0BnoDtGjRwoNIkRk+frKneLEofID7Lx/A/eC7502vHLvlC+pt3w50pMuWL1lcsznlTfG02/0zBZJDsFxjr7+TsQBd43Oilk66bPmS1mvXFSupVDoJSxZ1t2+HOlBje2TzZzIaG6P7+On2zB+qFzj32WDnljRLkny8KP1QTz1kU1VELsEx5RTtTOEprapOAaYA5OXlVQSrR8oRlJd7Xlt8/ko5deX60SkXpFsE3ymvdZ0IvVseRr35/0efbqn3dVQeGXL9YLZPGcs/uma+/yUvSn8lELgTcjOgzBwlETkJGAIcr1q8+/ZK4ISgtB/EI6hhpJLaNRtBmhp9ly55k81JNh/lHXUceUcdl7T8KyKjew8uPtYKMeIWH16U/hygrYi0An4DegEXB0YQkU7AY0B3VQ302PU2MFZE6rnnp4AzrdwwyjMH5x0GKd4ftogJV9lPpLyQiY7moip9Vc0Xkb44CjwXmKqqC0RkJDBXVWcAE4BawMuuo6LlqtpTVTeKyCicDwfAyKJBXSM2e2m6BnINw8gsPM3TV9WZwMygsGEBx2HdEqrqVGBqvAJmIrGob/Fj1xPDiJNkO+szUk/WO1wLps6OHQDU3F6xXS4bhh9kaw+z++EnUE83csaSn9Itiu9kpBuGRDh49Z/U2fJ/nFhr/3SLYhhGmuh42DH8ANA1/B4UFRVT+kHc0P+upORrneTMpdumT3lvj/KzatkwImFKP8WUB2dwhr88nwFb6BnZg9n000CL73/hkB3f0/7bH6LGzeT5woZhpB5r6acYFRgy+mGGgOsv3zAMI3VYSz9FZOccCKPCYy9uxmFKP0WYkcaoyNhYVOZgSj/FxPvjydb50ummlma+10UjuzCbfjnHVH36uO6716i66S/oemy6RTEM3zClbxhhuOumkekWwTB8x8w7hmGExUz5mYcp/RSx5/o/AWi23pyMGhURU/+ZQkYp/VM3fkyu5qdbjJCMvHEE/T56iruvGRJXevvJGYbhBxll03/6vL7pFiEig4c/EHsi0/aGYfhIRrX0DcMwjMh4Uvoi0l1EfhSRxSIyMMT140TkaxHJF5Hzg64ViMg892+GX4IbhmEYsRNV6YtILjAJOA1oD1wkIu2Doi0HLgdeCJHFdlXt6P71TFBewzBSyP5LVwGw98oNaZbE8AsvNv0jgMWqugRARKYBZwELiyKo6lL3WmESZMxqRJ0qrVxOB6iN8klt3cwWqZNwPiNuGskIgBM7JpyXUT7wYt5pCqwIOF/phnmlmojMFZHPReTsmKQzOKbtMZy24X/0+PyTdItiVCCu+fQN+nzzarrFMMohXlr6oTwBxDKnpIWqrhKR1sD7IjJfVX8pVYBIb6A3QIsWLWLIOvPpdtYZdEu3EEaF446h96RbBKOc4qWlvxJoHnDeDFjltQBVXeX+XwJ8AHQKEWeKquapal6jRo28Zm0YhmHEiBelPwdoKyKtRKQK0AvwNAtHROqJSFX3uCHQmYCxAMMwDCO1RFX6qpoP9AXeBn4ApqvqAhEZKSI9AUTkcBFZCVwAPCYiC9zkBwBzReRbYDYwXlVN6RuGYaQJUS1fSz7z8vJ07ty56RbDMAyjQiEiX6lqXrR4tiLXMAwjizClbxiGkUWY0jcMw8giTOkbhmFkEeVuIFdE1gHLEsiiIbDeJ3H8xOSKDZMrNkyu2MhEufZR1agLncqd0k8UEZnrZQQ71ZhcsWFyxYbJFRvZLJeZdwzDMLIIU/qGYRhZRCYq/SnpFiAMJldsmFyxYXLFRtbKlXE2fcMwDCM8mdjSNwzDMMJgSt8wDCOLyBilH23z9iSU11xEZovIDyKyQERucsPvEpHfAjaD7xGQZpAr348icmqyZBeRpSIy3y1/rhtWX0RmicjP7v96briIyES37O9E5NCAfC5z4/8sIpclKNN+AXUyT0Q2i8jN6agvEZkqIr+LyPcBYb7Vj4gc5tb/YjdtqI2IvMo1QUQWuWX/W0T2cMNbisj2gHqbHK38cPcYp1y+PTdx3LZ/4cr1kjgu3OOV66UAmZaKyLw01Fc43ZD2dwwAVa3wf0Au8AvQGqgCfAu0T3KZewOHuse1gZ9wNo6/C7gtRPz2rlxVgVauvLnJkB1YCjQMCrsHGOgeDwTudo97AG/i7JB2FPCFG14fWOL+r+ce1/Pxea0B9klHfQHHAYcC3yejfoAvgaPdNG8CpyUg1ylAJff47gC5WgbGC8onZPnh7jFOuXx7bsB0oJd7PBm4Pl65gq7fCwxLQ32F0w1pf8dUNWNa+sWbt6vqLqBo8/akoaqrVfVr93gLzl4DkfYOPguYpqo7VfVXYLErd6pkPwt42j1+Gjg7IPwZdfgc2ENE9gZOBWap6kZV/QOYBXT3SZZuwC+qGmnlddLqS1U/AjaGKC/h+nGv1VHVz9T5dT4TkFfMcqnqO+rsaQHwOc7OdWGJUn64e4xZrgjE9NzcFmpX4BU/5XLz/RvwYqQ8klRf4XRD2t8xyBzzTqKbtyeEiLTE2QbyCzeor9tNmxrQJQwnYzJkV+AdEflKnP2HAfZS1dXgvJTAnmmQq4helP4xpru+wL/6aeoe+y0fwJU4rboiWonINyLyoYgcGyBvuPLD3WO8+PHcGgCbAj5sftXXsbyOunYAAAKRSURBVMBaVf05ICzl9RWkG8rFO5YpSj/RzdvjL1ikFvAqcLOqbgYeBdoAHYHVOF3MSDImQ/bOqnoocBpwg4gcFyFuKuXCtdf2BF52g8pDfUUiVjmSVW9DgHzgeTdoNdBCVTsBtwAviEidZJUfAr+eW7LkvYjSDYuU11cI3RA2ahgZklJnmaL0E9q8PV5EpDLOQ31eVV8DUNW1qlqgqoXA4zjd2kgy+i67lmxG/zvwb1eGtW63sKhL+3uq5XI5DfhaVde6Mqa9vlz8qp+VlDbBJCyfO4B3BvB3tzuPaz7Z4B5/hWMvbxel/HD3GDM+Prf1OOaMSiHkjQs3r3OBlwLkTWl9hdINEfJL7Tvm1fhfnv+ASjiDHK0oGSQ6MMllCo4t7YGg8L0Djvvj2DcBDqT0ANcSnMEtX2UHagK1A44/xbHFT6D0INI97vHplB5E+lJLBpF+xRlAquce1/eh3qYBV6S7vgga2POzfoA5btyiQbYeCcjVHVgINAqK1wjIdY9bA79FKz/cPcYpl2/PDafXFziQ2ydeuQLq7MN01RfhdUP5eMcS/RGXlz+cEfCfcL7gQ1JQXhecLtV3wDz3rwfwLDDfDZ8R9OMY4sr3IwGj7X7K7r7Q37p/C4ryw7Gdvgf87P4venkEmOSWPR/IC8jrSpyBuMUEKOoEZKsBbADqBoSlvL5wuv2rgd04raar/KwfIA/43k3zMO7K9zjlWoxj1y16xya7cc9zn++3wNfAmdHKD3ePccrl23Nz39kv3Xt9Gagar1xu+FPAdUFxU1lf4XRD2t8xVTU3DIZhGNlEptj0DcMwDA+Y0jcMw8giTOkbhmFkEab0DcMwsghT+oZhGFmEKX3DMIwswpS+YRhGFvH/KE86Iw/LdqoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "continue_train = True  # set to True to continue to train\n",
    "num_big_steps = 50     # number of small steps\n",
    "num_small_steps = 200  # number of big steps\n",
    "if continue_train:\n",
    "    for k in range(num_big_steps):\n",
    "        for j in range(num_small_steps):\n",
    "            mine_dXY_list = np.append(mine_dXY_list, np.zeros((rep, 1)), axis=1)\n",
    "            mi_list = np.append(mi_list, np.zeros((rep, 1)), axis=1)\n",
    "            for i in range(rep):\n",
    "                minee_mine_list[i].step_mine()\n",
    "                mine_dXY_list[i,-1] = minee_mine_list[i].forward_mine()\n",
    "                mi_list[i,-1] = mine_dXY_list[i, -1].copy()\n",
    "        # To show intermediate works\n",
    "        for i in range(rep):\n",
    "            plt.plot(mine_dXY_list[i, :],label='dXY')\n",
    "            plt.title('Plots of divergence estimates')\n",
    "        display.clear_output(wait=True)\n",
    "        display.display(plt.gcf())\n",
    "    display.clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Save current results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current results saved.\n"
     ]
    }
   ],
   "source": [
    "overwrite = False  # set to True to overwrite previously stored results\n",
    "if overwrite or not os.path.exists(chkpt_name):\n",
    "    minee_mine_state_list = [minee_mine_list[i].state_dict() for i in range(rep)]\n",
    "    torch.save({\n",
    "        'dXY_list': dXY_list,\n",
    "        'dX_list': dX_list,\n",
    "        'dY_list': dY_list,\n",
    "        'mine_dXY_list': mine_dXY_list,\n",
    "        'mi_list': mi_list,\n",
    "        'minee_mine_state_list': minee_mine_state_list\n",
    "    }, chkpt_name)\n",
    "    print('Current results saved.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the ground truth mutual information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ground truth is 0.4084425498386879 nats.\n"
     ]
    }
   ],
   "source": [
    "mi = mg.ground_truth * d\n",
    "print('Ground truth is {} nats.'.format(mi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply moving average to smooth out the mutual information estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "mi_ma_rate = 0.01            # rate of moving average\n",
    "smooth_mi_list = mi_list\n",
    "for i in range(1,smooth_mi_list.shape[1]):\n",
    "    smooth_mi_list[:,i] = (1-mi_ma_rate) * smooth_mi_list[:,i-1] + mi_ma_rate * smooth_mi_list[:,i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the mutual information estimate after different number of iterations. The red dashed line shows the ground truth, and the green dotted line is the number of iterations where 90% of the ground truth is reached."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAELCAYAAADkyZC4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4FdXWwOHfItK7gIoUAUUREHIhgFgARYrUq6JgQVCv6EWuehU/Gzawd71WVBARBRUVkCaCYAGEIIgUkaqE3lsKKev7YyaHk+QkOQmZMyFZ7/PkYcqemZVJOCuz9569RVUxxhhjCloJvwMwxhhTNFmCMcYY4wlLMMYYYzxhCcYYY4wnLMEYY4zxhCUYY4wxnvA0wYhIVxFZIyLrROSBHMr1EREVkRh3vZ6IJIjIMvfrHS/jNMYYU/BO8urEIhIFvAl0AuKAxSIyWVVXZSpXEbgT+CXTKdararRX8RljjPGWl08wrYF1qrpBVY8C44HeIcqNAJ4HEj2MxRhjTIR5mWBqAZuD1uPcbQEi8g+gjqp+E+L4+iKyVETmicjFHsZpjDHGA55VkQESYltgXBoRKQG8AgwMUW4bUFdV94hIS+BrEWmiqgczXEBkEDAIoHz58i0bNWpUULEbY0yxsGTJkt2qWsOLc3uZYOKAOkHrtYGtQesVgabAXBEBOA2YLCK9VDUWSAJQ1SUish44G4gNvoCqjgRGAsTExGhsbIbdxhQqmw84D/R1KtfJpaQxkSMif3l1bi8TzGKgoYjUB7YA/YDr0neq6gGgevq6iMwFhqpqrIjUAPaqaqqINAAaAhs8jNUYz/X/qj8AcwfO9TcQYyLEswSjqikiMgSYCUQBo1R1pYgMB2JVdXIOh7cDhotICpAK3K6qe72K1ZhIGNZumN8hGBNRUlSG67cqMmOMyTsRWaKqMV6c297kNyZCNuzbwIZ9VtNrig8v22CMMUFunnQzYG0wpviwBGNMhDzR4Qm/QzAmoizBGBMh7eu19zsEYyLK2mCMiZA1u9ewZvcav8MwJmLsCcaYCLntm9sAa4MxxYclGGMi5OmOT/sdgjERZQnGmAi5oM4FfodgTERZG4wxEbJi5wpW7FzhdxjGRIw9wRgTIUOmDQGsDcYUH5ZgjImQFzq94HcIxkSUJRhjIqRVrVZ+h2BMRFkbjDERsmz7MpZtX+Z3GMZEjD3BGBMhd8+4G7A2GFN8FJ0Es2YNdOiQcds118DgwRAfD926ZT1m4EDna/du6NMn6/5//xv69oXNm6F//6z7770XevZ0rn3bbVn3DxsGl10Gy5bB3Xdn3f/003DBBTB/Pjz0UNb9r74K0dHw3Xfw5JNZ97/7LpxzDkyZAi+9lHX/2LFQpw5MmABvv511/xdfQPXq8OGHzldm06ZBuXLw1lvw2WdZ98+d6/z74ovwzTcZ95UtC9OnO8sjRsDs2Rn3V6sGEyc6yw8+CAsWZNxfuzZ8/LGzfPfdzj0MdvbZMHKkszxoEPz5Z8b90dHO/QO44QaIi8u4v21beOYZZ/mqq2DPnoz7O3aERx5xli+/HBISMu7v0QOGDnWWM//eQcjfvVcrHHb2fdjBfvfsd89ZjtDvXga5/e4VIE+ryESkq4isEZF1IvJADuX6iIiKSEzQtgfd49aISBcv4zQmEqIPVyD6cAW/wzAmYjybcExEooA/gU5AHM4Uyteq6qpM5SoCU4FSwBB3yuTGwKdAa+B04DvgbFVNze56NuGYKewWb1kMWGO/KVxO1AnHWgPrVHWDqh4FxgO9Q5QbATwPJAZt6w2MV9UkVd0IrHPPZ8wJ675Z93HfrPv8DsOYiPGyDaYWsDloPQ5oE1xARP4B1FHVb0RkaKZjF2Y6tpZXgRoTCW90e8PvEIyJKC8TjITYFqiPE5ESwCvAwLweG3SOQcAggLp16+YrSGMipekpTf0OwZiI8rKKLA6oE7ReG9gatF4RaArMFZFNwPnAZLehP7djAVDVkaoao6oxNWrUKODwjSlY8zfPZ/7m+X6HYUzEePkEsxhoKCL1gS1AP+C69J2qegConr4uInOBoW4jfwLwiYi8jNPI3xBY5GGsxnjuodlOd2B7D8YUF54lGFVNEZEhwEwgChilqitFZDgQq6qTczh2pYh8BqwCUoA7cupBZsyJ4N0e7/odgjlOqWlpHElMoVK5Un6HckLwrJtypFk3ZWOM157/ehmzf9/ChHsuo0r50n6HE1JyahoAJaPCawE5UbspG2OCzNs0j3mb5vkdRqF1JDGZtAL6g3fz7sMkJRd8pcfs37cA0Pfl71BVuoyYysxlxzrLpqny+YL1HE0Jfe2D8UdJcRNAKDsPJPD2zJU5lsnJnkOJ9Hh6Oj2enh5y/4H4o3ww+w+mLvkrX+fPK0swxkTIY3Mf47G5j/kdRqGSnJrGwfijHElK5soXvuWWt+Ye1/nmrdzK69N+519vz6PXszMCH/RHkpIZOWsV+w4nZfnwTk1zEsWXv2wMec4uI6bS7+XvsmzfccAZwuXlKcvpMmIqAB/M/oP3v/uDns/MCJRbsn4Xz321lCUbdnH1S7Po/vT0bBNQ/9fn8PWiTTz0Segm5/lrtrPzQELIfQDXvXpsWJz073P/kSTu+2gBicmpXPPSLD6bv57Xp61gx/74bM9TUIrOWGTGFHKjeo/yO4Q82bTzEE9+sYSXB15QIG0Oaapc/uQ02p59KnVrVODmSxtx0xvfs+tgIo1qVQFg6954+r8+h4/+cwkiod5WyN60X//mtam/Z9jW85kZzHykO1c+/y0AExc6SWTmI90DZbo9NQ2Ad79dRZM6ValSrhSnVikHQHoTwr4jSSz/69iYYV3/UYcB//s+w7VS05QvFmwIrHcZMZU3/nVRIFnMWXGsI2x6Apo+rBslRNh1MIE3p68M7P9t057AOQCuv7gh435cG9j/6k0X8NAnixh9R4dsq+q6Z3qK6f3sjAzr67YfDHyfXrEnGGMipEHVBjSo2sDvMHJ0IP4or3yznN827eG2d39g854jDP98SbblVZWc2nG7jJhK1xFTUVXe/241AAv+3MGEn9cDsOugM4DHH1v2B47ZeSCB0d+vITE5ldS0ND6Y/QddRkwNfEDe/u4PdBkxlcXrdtJlxFQOJhwFyJJc0oV6WkhNS6PLiKl8EvShDXDnBz9z4/++58+tTjwvTv4tsO++j469+z1j6WYy+3Lhhizbhrz/U8iY0t31wc90GTGVG16bw4I/d2TYN3rOH4HlcZnivHv0fOKTUrju1dn8vesQU5f8FUhG4dq+Pz7Px+SVNfIbEyHfbXCqWS5rcFlEr/vKN8uZsXQzXaJrc0/P5oHt81ZuZdbyOBav28XUhy5n/Y6D3PnBzyHPcUWb+vSMOYOD8Uf5bnkc5599KsM+XRzY/8ndHalWsUyGY3YeSKD/63MC6xc2Oo2f/9geWH/smpY88Vn2yQsg5swaxK7fFVjv1eoMJi/O2n7w0X8u4cZMTxTpTqlcNsdqpezMfKQ7t73zA5t2HcrzsX7KfM9y8+2jPTxr5LcEY0yEdPiwA5D9ezCqmudqoVCOpqTS85kZvHrTBZxbu6rnf6UC3NOzGdv2xVO1Qml6t6oHwNptB3L9C94LT17bKkPyy6+WDaqzZMPuAogoo/t6N+eFSb/lXjCfXh7Ylns+XJDt/ru6n5fhac/LBGNtMMZEyNgrxobc/sj4xSxauxOADwa3p3a10EP6p6Yp3Z5y2jDu7dWcimVLArB++0EOJybTvF41AJ6euBRwqlHu/2d0QX8bIb08ZXlg+a0ZK3Mo6b1WZ51CvwvPZLxbDZfu6rYN+HxB1mqs7GSXXLq1qMu0X//O8diZj3Rn/faDDH7vxyzbATqeV4uuT04LO5Z0t3VuzLvfrsqxTM2q2berDLzkHLpE18m2OrGgWYIxJkLqVK6TYT01TdlzKDGQXABueWse/+rYiPdnO/Xv6R9IR5KSAw3VC/7cQZ8Xv2XgJeewccdB5q3alqFscF3+c1/nf4rmKQ92zdAbqqDUOrk8W/YeCax/MbQzqsqidTuP+y/7z4d2AuCmSxuRcDSVSYs3Bfb967JzqVi2FKOC2jbC9Vz/Nvy16zC9W9Vj5rLNgQRze+fGXNGmfoanxHcGXQzAmadVCjxN3NOzGe0a1wyUERGmPNiVPYeSGPjGsaq9ZmeczFPXteb2d38M3KMXbjyf+z5aSP92DbmyTf1Aghl3V0eufy3TZGrAyRXKMO3hbsTtOcwZNSqy/K89nHlaJcqXLpmlbLcWdfk2z3cjfFZFZkyEzFjnfFh3Pasr+48k0TdE19fMZj7Snd//2sPQjxbmWrYgjOjXiokLN3DNhWfSskENFqzZweOfHd//q7duvYjB7zlVZekvKKZ/IHdrUZe7up+XofzBhKNc/eKsDNua16sW6Fl1ds3K/LntAACNa1dlVdy+QLng3mHBSXnGsG6ICKoaeHKoUalMoJNB5XKleH9we0qdFMWmnQe5a9SxMeOe69+G6HqBUa0CbVrB1wtOMMExQO5VnxN+Xs+oOX/w5f91DiSB5NQ0ejw9nWFXteDixjVZ8fdezq1dhagSGftl/bJ2ByVEiK5fnVm/xbH7YCI3djg722ul+3Prfpb/tZc+bRt4+qKlJRhjIiS9DWZk56/5TzaN6V65pWMjPpgd+i/3qBJCaprzOZD5wxHg+xVbePYr50noi6GdeXvmSmb/voVpD3cLdPHNybSHLyc5JY0SJYRSJ0UBxz6QQ10veP/AS86hTcNTqFi2JDe8Nidwvvs+WsjKzfsytCeE6mjwyjfLufwfdWhUq2pg2/rtBxj83k+8e1s7Tq5Ymue+WsawPi0oW+pYhc6TXyzhx9VOh4QpD3YNxA0Qn5TCFc/PpO8FZ3Jzx0bAsYSQ0/dUWFmCCYMlGFPYdRjuzPNeWqrmUrLgzRjWLfCX+3u3t6NECWH/kaPcO2YBT13Xmka1qrDnUCJn1KgY8vicEkJOnQiy+7BVVRQokc1f9uN/Wsfo79cEjk9/hyb9nGmqLFq7kzYNT0FEOJqSmiEJ+OHHVdtISUvjkqYn1tRVlmDCYAnGFHa59ebK3DaRnR4t6/LNkmONzHWqlWfznpyPm/lId/7adQgRoW71Y50IjiQlh6ybz2zPoURKl4yiQpmsZWPX7+LhTxbxRN8Yzj/7VNZtO8Db366id6t6Gdodjlf6ezHhxGvCZwkmDJZgTGE27de/efib9wE4RZyJXW++9BxGzVnDPT2b0SXa6QDw9MRfmbdqW7YN7IO7NKZ7yzOY/fsWLmtWm6gSzhPAqrh9/Hd0xrlmLm16Oj1b1eOc0ytnqbs3Jp2XCcZ6kRkTAZt2HuIv/QpwEszrt1zIOadXoe+FZ2Uo99BVLXjoKmf5gnNOZf6aYz3CSp9Ugl6t6iEigYSUrnHtqrx160XUqlaBg/FHqVK+lO9VRsZYgjEmAqLrV6P5ogd59JqWXNLonLCOeeyamEC12ksD2tK07sk5lj/ztMoAlKlc9viCNaaAWIIxxmMzl21m3I9rKSWVaVgjbw3AJQTSlFyTizGFkacVsyLSVUTWiMg6EXkgxP7bReR3EVkmIj+JSGN3ez0RSXC3LxORd7yM0xgvvTxlOTv2J7BD5zNjfbYTuYb0+dDOjLuro0eRGeMtz55gRCQKeBPoBMQBi0VksqoGj3Pwiaq+45bvBbwMdHX3rVfVyIxzYYwHvl60kYPxyYH1v3UyH6/6mX+1vi7sc1QoUzJkzy1jTgReVpG1Btap6gYAERkP9AYCCUZVDwaVLw8UjS5txgBvz8w4ZtRVNV/gxX5tfYrGmMjzMsHUAoInTYgD2mQuJCJ3APcApYBLg3bVF5GlwEFgmKr+GOLYQcAggLp16xZc5Mbk0+bdh5m/Zjubd2d9L6V6+apULlPZh6iM8YeXCSbUK7pZnlBU9U3gTRG5DhgGDAC2AXVVdY+ItAS+FpEmmZ54UNWRwEhw3oMp6G/AmNx8/MNa2jWuyYYdB3nmy6U5ll15cBYTVmyib9O+EYrOGH95mWDigODO+rWBrdmUBRgPvA2gqklAkru8RETWA2cD9ialKRRmL4/jeXfk37Hz/sy23PRh3fjftBVM+/Vvfto5gYOx1SzBmGLDywSzGGgoIvWBLUA/IEPrpog0VNX0uUC7A2vd7TWAvaqaKiINgIZA+BM5GFPAUtOU1LQ0Sp0UleuQLzUqlaFWtfIkHk2lhEhgaPcW8jiTr788EuEaUyh4lmBUNUVEhgAzgShglKquFJHhQKyqTgaGiMhlQDKwD6d6DKAdMFxEUoBU4HZV3etVrMbk5I8t+zIM3x6sbvUK/L37cGA91Ii+6U6tVIVyJbOfDMqYosbGIjMmhCmxf3FSlHAoITnbYe5fuPF8mp3hzCJ5JCmZ/UeOUuvk8lnKrdt2gDve/4l/dt7pDDvf7AZPYzcmL2wsMmMi4HBiMmVLncS/3p7L1r3xOZaddH8XygTNH1K+dMlsR/k9q2ZlZj7SPTAfjCUYU1xYgjHF3v1jF7LMnS0xO1e3bUDl8qW44JzTQj6lhGNW/1m5FzKmCLEEY4qlg/FHmbdqK29MX5lr2ceviaHtOace9zVLRtkb+aZ4sQRjigVV5fe/9/LL2p18sSD3DolTHuzK73/tpUmdqhmqwo7Hh8s+BGBg9MACOZ8xhZ0lGFPkqWpguuDsdGpem8FdmlCu9LH/Ei3PrFGgcViCMcWNJRhTJKUnlbKlokg4mhqyTN8Lz2TCz+uZ8mDXiEzONXfgXM+vYUxhYgnGFEnpTyyhksvgLo3p3bo+ADdf2iiicRlTnFiCMSe0NFVKiKCqvDF9BalpyvSlm0OWnfxAV0qX9G8a4feWvAfArS1v9S0GYyLJEow5YfV6ZjpJKWnUOrk8FcuW5I8t+7OU+fiuS0lOSeP0fHYtLkgTVk4ALMGY4sMSjDkhrdy8l6SUNAC27M06ND74/8SS2Xc3fud3CMZElCUYc0LIbYDJYM/d0Ibo+tU9jMYYEw5LMKbQmxK7KaxyHwxuz8GEZBrXruptQPn01uK3ABjcarDPkRgTGZZgTKH2zJdLmbsy+2mEZj7SPYLRHJ8pf04BLMGY4sMSjCkUnKFbttGjZV1EJGSV2Lm1qrB6y36anXEyy//ay2NXt/Qh0vybfv10v0MwJqIswRhf/bXrEIPe+SGw/sb0FSHLPXBFNB2anI5IqJm4jTGFUQkvTy4iXUVkjYisE5EHQuy/XUR+F5FlIvKTiDQO2vege9waEeniZZzGH4nJqRmSS3ZG3dGBS5rWOuGTy2sLX+O1ha/5HYYxEePZE4yIRAFvAp2AOGCxiExW1VVBxT5R1Xfc8r2Al4GubqLpBzQBTge+E5GzVTX0mB/mhKGqiAiJR1Po/dzMHMu+dvMFNKpVOBvs82P2xtkA3HX+XT5HYkxkeFlF1hpYp6obAERkPNAbCCQYVT0YVL48kD69Zm9gvKomARtFZJ17vgUexms8kpicyra9R7h95I8h9591WiXWbT/IZ/d24tOf1rE6bh+DOp1bpJILwORrJ/sdgjER5WWCqQUEj9kRB7TJXEhE7gDuAUoBlwYduzDTsbW8CdN4rfezM7Ldd/G5NRnWp0Vg/fbOjbMta4w5sXjZBhOqwlyzbFB9U1XPBO4HhuXlWBEZJCKxIhK7a9eu4wrWFLyjKam5viAZnFyKuhfnv8iL81/0OwxjIsbLBBMH1Alarw1k/0IDjAf+mZdjVXWkqsaoakyNGgU7d4c5PqlpSs9nMj65lC4ZxcxHunPWaZUAZ5yw4mRB3AIWxFktryk+RDXLg0HBnFjkJOBPoCOwBVgMXKeqK4PKNFTVte5yT+AxVY0RkSbAJzjtLqcDs4GGOTXyx8TEaGxsrCffiwlfmioL1uxg+OdLMmx/7/Z21K1R0aeojDHZEZElqhrjxbk9a4NR1RQRGQLMBKKAUaq6UkSGA7GqOhkYIiKXAcnAPmCAe+xKEfkMp0NACnCH9SAr/J74LJb5a3Zk2FZCYPqwE+dte2NMwfHsCSbS7AnGX8v/2sN9Hy3Msn3GsG4n/PsrBeXZn54F4IGLsrwSZoxvTsgnGFP07TmUyA+rtjFm7posM0fWqVae9wd38CewQmrZ9mV+h2BMRFmCMfmSpsp1r84Oue/d29pR7xRrb8lsfJ/xfodgTERZgjF5tu9wEv1eyTp51qs3XcC5hXSofGNM5FmCMXmy/0jW5FK5XCk+u7eTTxGdOEbMGwHAI+0f8TkSYyIj1wQjImcDbwOnqmpTEWkG9FLVJz2PzhQqhxKS6ftyxuRyIs3H4rc1e9b4HYIxERXOE8x7wH3AuwCqulxEPgEswRQjOw8k0P/1ORm2tW54ik/RnJg+vvJjv0MwJqLCSTDlVHVRpq6mKR7FYwqhUMO9vPGviwJv5BtjTCjhJJjdInIm7lhgItIH2OZpVKbQGPbpogzrt3duzBVt6vsUzYnt0e8fBWD4JcN9jsSYyAgnwdwBjAQaicgWYCNwvadRmUJj8bpjg4he2vR0Sy7HYfPBzbkXMqYICSfBqKpeJiLlgRKqekhE7FOmiMs87Iu9kX/8Rvce7XcIxkRUOKMpTwRQ1SOqesjd9oV3IRm/qWqG5PLwVS0suRhj8izbJxgRaYQzZXFlEbkyaFcloIzXgRn/3Dnq5wzr7RrX9CmSouXB7x4E4JnLnvE5EmMiI6cqsnOAHkAVoGfQ9kPArV4GZfwT3GPs3l7N6Ny8Tg6lTV7sSdjjdwjGRFS2CUZVJwGTRKStqtosScVA5u7IllwK1sieI/0OwZiICqeRf6mI3IFTXRaoGlPVmz2LykRMUnIqr36znDkrjk0Yen7DU3iiXysfozLGFAXhNPKPBU4DugDzcKYvPpTjES4R6Soia0RknYhkmQRDRO4RkVUislxEZovIGUH7UkVkmfs1Obxvx+RV/9fnZEgugCUXjwz9dihDvx3qdxjGREw4CeYsVX0EOKKqY4DuwHm5HSQiUcCbwOVAY+BaEWmcqdhSIEZVm+H0THs+aF+Cqka7X73CiNPk0bvfruJA/NEM26Y93M2naIq+hOQEEpIT/A7DmIgJp4os2f13v4g0BbYD9cI4rjWwTlU3AIjIeKA3zjTIAKjq90HlFwI3hHFeUwASj6bw5S8bA+s2aKX33uz+pt8hGBNR4SSYkSJSFXgEmAxUAB4N47haQPCry3FAmxzK3wJMD1ovIyKxOOOePauqX4dxTZOLu0f9zOot+zNss+RijPFCrglGVd93F+cBDfJw7lBv5mnIgiI3ADFA+6DNdVV1q4g0AOaIyO+quj7TcYOAQQB169bNQ2jF0/0fL8ySXCbcc5lP0RQ/d8+4G4BXu77qcyTGREY488FUAW7EqRYLlFfVO3M5NA4I7udaG9iauZCIXAY8DLRX1aSg8291/90gInOBfwAZEoyqjsQZJ42YmJiQycs4+r48i/1HMra3VK9UhirlS/sUkTGmqAunimwaTvvI70BaHs69GGjojlu2BegHXBdcQET+gTPPTFdV3Rm0vSoQr6pJIlIduJCMHQBMHgUnl0/u7kgJEapWsOQSSfbkYoqbcBJMGVW9J68nVtUUERkCzASigFGqulJEhgOxqjoZeAGnTedzd6yrv90eY+cC74pIGk5Pt2dVdVXIC5lcBb9AOX1YN0rYuGLGmAgIJ8GMFZFbgW+A4CqsvbkdqKrTcJ6Agrc9GrQcsgFAVecTRldok7P/fPATf249kGGbJRf/3DH1DsB6k5niI5wEcxTnSeNhjjXSK3lr8DcRdjgxOUty+er/uvgUjQEoW7Ks3yEYE1HhJJh7cF623O11MKZgZG7Qj65fjeduON/HiAzAi51f9DsEYyIqnASzEoj3OhBTMFbH7cuQXD6+61JqVLK/nI0xkRdOgkkFlonI92Rsg8mtm7KJsKTkVO4ePT+wPuTyJpZcCpFBUwYBNqqyKT7CSTBfu1+mENu08xC3vftDYH3aw5cTVSKcoeZMpFQrW83vEIyJqHDe5B8TiUBM/k2J/Ys3pq8IrA/t1dySSyFkM1ma4ianKZM/U9VrROR3Qgzx4o6AbHz24+ptGZJLzJk16NS8to8RGWOMI6cnmLvcf3tEIhCTP09+8Wtg+YwaFXjqutY+RmNyctOkmwAY3Xu0z5EYExk5TZm8zV0crKr3B+8TkeeA+7MeZSJFVen65LF3WJ+8thUxZ9bwMSKTmzqVbApqU7yE08jfiazJ5PIQ20wEBScXgFZnneJTJCZcwy8Z7ncIxkRUTm0w/wYGA2eKyPKgXRWBn70OzGTv3jELAstdo+twZ/emPkZjjDGh5fQE8wnOBGDPAA8EbT8UzjhkpuClpindnsr45PLfntbX4kRxw5fOhK0fX/mxz5EYExk5tcEcAA6IyDBguzt0fgegmYh8pKr7szvWFLzgEZHTzRjWzYdITH6dU+0cv0MwJqLCaYOZCMSIyFnABzjTJn8C2KdbhNw/dmGG9ZpVy/HhkEt8isbk1yPtH/E7BGMiKpwEk+bO7XIl8Kqq/k9ElnodmDlm2aY9geXrLj6LAR3sL2FjTOEXzuveySJyLc60yd+420qGc3IR6Soia0RknYg8EGL/PSKySkSWi8hsETkjaN8AEVnrfg0I53pF0ZNfLAksz3ykuyWXE1i/L/rR74t+fodhTMSE8wRzE3A78JSqbnSnQM61lVJEooA3cbo5xwGLRWRyppkplwIxqhrv9lp7HugrIicDjwExOKMILHGP3ZeXb+5E98vaHfy4ejsA/S480+dozPGKPi3a7xCMiahwxiJbJSL3A3Xd9Y3As2GcuzWwTlU3AIjIeKA3EEgwqvp9UPmFwA3uchdgVnpvNRGZBXQFPg3jukXGo+NjA8s3XdrIx0hMQXjgoiwP8cYUablWkYlIT2AZMMNdjxaRyWGcuxawOWg9zt2WnVtwukXn59gi56O5fwaWZz7S3cdIjDEmf8Jpg3kc52lkP4CqLgPqh3FcqMnfswyaCSAiN+BUh72Ql2NFZJCIxIpI7K5du8II6cSQeDSFcT+uBaBCmXBqMc2J4KrPruKqz67yOwxjIiacBJPivhMTLGSiyCQOCB58qTawNXMhEbkMeBjopapJeTmZ2Cs1AAAgAElEQVRWVUeqaoyqxtSoUXTG4Zowf31g+YuhnX2MxBSktrXb0rZ2W7/DMCZiwvnzeIWIXAdEiUhD4E5gfi7HACwGGrqdArYA/YDrgguIyD+Ad4GuqrozaNdM4GkRqequdwYeDOOaJzRV5ekvl/LDKmec0Qn3XIZIqIc5cyIaesFQv0MwJqLCSTD/wXnCSMJ5wXIm8GRuB7nvzgxxy0cBo1R1pYgMB2JVdTJOlVgF4HP3g/RvVe2lqntFZAROkgIYXhyGp8k8gGWV8qV9isQYY46fqIZT21X4xcTEaGxsbO4FC7Hg4WDG3nkpp1Qu62M0pqD1+rQXAJOvDaePjDGRISJLVDXGi3NbC3Ih8MeW/dw16tgA1ff0bGbJpQjqWL+j3yEYE1GWYHx2/Wuz2X0wMbD+z9b16BJtE1MVRXedf1fuhYwpQsLpRWY8oqoZkgvAv7s08SkaY4wpWDlNOPY/cuiOrKp3ehJRMfLS5GPzuM0Y1s16jBVxl4+7HIDp10/PpaQxRUNOVWQndot5Ibc6bh+zlscB8PR1rS25FAM9z+7pdwjGRFROE46NiWQgxUmaKnePPvYqUcszi85LoiZ7g1sN9jsEYyIqpyqyHPtSqmqvgg+neLg86H2Xz+/t5GMkxhjjnZyqyNriDDj5KfALoccHM3mUeerjSuVK+RSJibTLProMgO9u/M7nSIyJjJwSzGk4c7lcizPEy1TgU1VdGYnAiqLg5HJyhdK8eKONS1Wc9G3S1+8QjImonNpgUnGG6J8hIqVxEs1cERmuqv+LVIBFQXxSCre8NTewfuaplXhr0MX+BWR8cWvLW/0OwZiIyvFFSzexdMdJLvWA14EvvQ+r6EhNS+OK52cG1v/dpTH/bB3ObAfGGHNiy6mRfwzQFGcSsCdUdUXEoipCgt91ASy5FGMdPuwAwNyBc32Nw5hIyekJpj9wBDgbuDPoPQ0BVFUreRxbkTD79y0AvHtbO+qdUtHnaIyfBkYP9DsEYyIqpzYYG0bmOP26YXdg2ZKLsQRjihsb7NIj36/YwrNfLQPg6etb+xyNKQySU5MBKBlV0udIjIkMT59SRKSriKwRkXUi8kCI/e1E5FcRSRGRPpn2pYrIMvfrhJpAo8uIqYHkAtCygb2pb6DT2E50Gmsv1priw7MnGBGJAt7EeZcmDlgsIpNVdVVQsb+BgUCouWQTVDXaq/i8sn77gQzrnw+1DxTj+FeLf/kdgjER5WUVWWtgnapuABCR8UBvIJBgVHWTuy/NwzgiZsn6XTz0ySIAqpQvxYR7LLmYY25odoPfIRgTUV5WkdXCGWomXZy7LVxlRCRWRBaKyD8LNrSCt377wUByARj/38t8jMYURvHJ8cQnx/sdhjER4+UTTKixy7KdXyaEuqq6VUQaAHNE5HdVXZ/hAiKDgEEAdevWzX+kxyk5NY3B7/0YWH/xxvNt+H2TRbdx3QB7D8YUH14mmDggeO7f2sDWcA9W1a3uvxtEZC7wD2B9pjIjgZEAMTExeUleBarH08cmkGp79qmcd0Y1v0Ixhdi/Y/7tdwjGRJSXCWYx0FBE6gNbgH44g2bmSkSqAvGqmiQi1YELgec9i/Q4BA9gOe3hbkSVsCcXE1rfpjbYpSlePGuDUdUUYAgwE1gNfKaqK0VkuIj0AhCRViISB1wNvCsi6SM1nwvEishvwPfAs5l6nxUKi9buDCy3a1zTkovJ0YHEAxxIPJB7QWOKCFH1rWapQMXExGhsbORmeT6UkEyfF78FoFrF0nxytzXqm5zZWGSmMBKRJaoa48W57U3+fEhTDSQXgP/dcpGP0ZgTxZ1t7vQ7BGMiyhJMHqWmpdHtqWON+l/f34Wypew2mtxdee6VfodgTETZgJZ5FJxcRt7ezpKLCdvu+N3sjt+de0Fjigj7dMyDlNRjAw7c3rkxZ9SwEZJN+Pp85gy3Z20wpriwBJMH3YPed7mijU0cZvLm3rb3+h2CMRFlCSZMi9cd65L8RF9POlyYIq7nOT39DsGYiLI2mDBs2nmIYZ8uBuDCc07l/LNP9TkicyLafng72w9v9zsMYyLGnmBycTD+KLe9+0Ng/dFr7OnF5E+/L/oB1gZjig9LMDlITVOufmlWYP3ze234fZN/D1yUZc49Y4o0SzA5+GXtjsDylAe7UuqkKB+jMSe6rmd19TsEYyLK2mBy8NUvGwFnbhdLLuZ4bT6wmc0HNude0Jgiwp5gcrD8r70AVK1Q2udITFHQ/6v+gLXBmOLDEkw2lm60N65NwRrWbpjfIRgTUZZgsvHDqm0AvHrTBT5HYoqKyxrYiNumeLE2mGys2ryP+qdU5NzaVf0OxRQRG/ZtYMO+DX6HYUzEWIIJYe6KrWzadYjLmtX2OxRThNw86WZunnSz32EYEzGeJhgR6Soia0RknYhkeQlARNqJyK8ikiIifTLtGyAia92vAV7Gmdn8Nc7b1s3rVYvkZU0R90SHJ3iiwxN+h2FMxHjWBiMiUcCbQCcgDlgsIpMzTX38NzAQGJrp2JOBx4AYQIEl7rH7vIo32Dy3/eWs0ypF4nKmmGhfr73fIRgTUV4+wbQG1qnqBlU9CowHegcXUNVNqrocSMt0bBdglqrudZPKLCDib6mJSKQvaYqwNbvXsGb3Gr/DMCZivOxFVgsIfqssDmhzHMfWylxIRAYBgwDq1q2bvygzUdUCOY8xmd32zW2AvQdjig8vE0yoP//D/fQO61hVHQmMBIiJiSmQzHAoIRmALtHWwG8K1tMdn/Y7BGMiyssEEwfUCVqvDWzNw7EdMh07t0CiysWGnQcBiK5XPRKXM8XIBXXsnSpTvHjZBrMYaCgi9UWkFNAPmBzmsTOBziJSVUSqAp3dbZ7bsT8BgNOqlovE5UwxsmLnClbsXOF3GMZEjGdPMKqaIiJDcBJDFDBKVVeKyHAgVlUni0gr4CugKtBTRJ5Q1SaquldERuAkKYDhqrrXq1iD7T6YCMCZp1oPMlOwhkwbAlgbjCk+PB0qRlWnAdMybXs0aHkxTvVXqGNHAaO8jC+UbfvjqVaxNKVL2ujJpmC90OkFv0MwJqJsLLJMtu+L57QqVj1mCl6rWq38DsGYiLKhYjLZvt8SjPHGsu3LWLZ9md9hGBMx9gQTJDk1jd0HEy3BGE/cPeNuwNpgTPFhCSbIrgMJKHBa1bJ+h2KKoFe7vup3CMZElCWYIDsPOl2UT6lkCcYUvOjTov0OwZiIsgQTZN/hJABOtimSjQcWb3F63ReHxv7k5GTi4uJITEz0OxTjKlOmDLVr16ZkyZIRu6YlmCB73QRTtUIZnyMxRdF9s+4DikcbTFxcHBUrVqRevXo2aGwhoKrs2bOHuLg46tevH7HrWoIJsvdwEiWjSlChjN0WU/De6PaG3yFETGJioiWXQkREqFatGrt27Yrode2TNMiBI0epUr6U/acwnmh6SlO/Q4go+39UuPjx87D3YIIcTkymQpnI1U+a4mX+5vnM3zzf7zBMhMydO5cePXpk2b5s2TKmTZsW4ojcPf30sRG5N23aRNOmhfuPFkswQY4kJVPeEozxyEOzH+Kh2Q/5HYYJkpKSEvFr5pRgcosnOMGcCKyKLMiRxBRqVLIGfuONd3u863cIxcqIESMYN24cderUoXr16rRs2ZKhQ4fSoUMHLrjgAn7++Wd69epFnz59uPnmm9m1axc1atRg9OjR1K1bl4EDB9KjRw/69OkDQIUKFTh8+DBz587l8ccfp3r16qxYsYKWLVvy8ccfIyLMmDGDu+++m+rVq9OiRYssMR09epRHH32UhIQEfvrpJx588EFWr17N1q1b2bRpE9WrV6dz587ExsbyxhtOm12PHj0YOnQoM2bMICEhgejoaJo0acJTTz1Famoqt956K/Pnz6dWrVpMmjSJsmULz2sWlmCCHE5Kpl6Zin6HYYqoc6qf43cI/unQIeu2a66BwYMhPh66dcu6f+BA52v3bnA/5APmzs3xcrGxsUycOJGlS5eSkpJCixYtaNmyZWD//v37mTdvHgA9e/bkxhtvZMCAAYwaNYo777yTr7/+OsfzL126lJUrV3L66adz4YUX8vPPPxMTE8Ott97KnDlzOOuss+jbt2+W40qVKsXw4cMzJJDHH3+cJUuW8NNPP1G2bFk+/PDDkNd89tlneeONN1i2zBluaNOmTaxdu5ZPP/2U9957j2uuuYaJEydyww035Bh7JFkVWZAjicmUtx5kxiPzNs1j3qZ5fodRLPz000/07t2bsmXLUrFiRXr27Jlhf/CH/4IFC7juuusA6N+/Pz/99FOu52/dujW1a9emRIkSREdHs2nTJv744w/q169Pw4YNEZE8fdD36tUrX08e9evXJzraeYG3ZcuWbNq0Kc/n8JJ9mrrSVIlPSqFCaWuDMd54bO5jQPF4DyaLnJ44ypXLeX/16rk+sWSmmvMM6uXLl892X3pvq5NOOom0tLTA+Y4ePRooU7r0sZexo6KiAm0n+e2pFRxP8HWBHF9WzRxHQkJCvq7vFU+fYESkq4isEZF1IvJAiP2lRWSCu/8XEannbq8nIgkissz9esfLOAESjqaQplgjv/HMqN6jGNU74lMcFUsXXXQRU6ZMITExkcOHDzN16tRsy15wwQWMHz8egHHjxnHRRRcBUK9ePZYsWQLApEmTSE5OzvGajRo1YuPGjaxfvx6ATz/9NGS5ihUrcujQoWzPU69ePZYtW0ZaWhqbN29m0aJFgX0lS5bMNY7CxLMEIyJRwJvA5UBj4FoRaZyp2C3APlU9C3gFeC5o33pVjXa/bvcqznRHEp2/QKyKzHilQdUGNKjawO8wioVWrVrRq1cvmjdvzpVXXklMTAyVK1cOWfb1119n9OjRNGvWjLFjx/Laa68BcOuttzJv3jxat27NL7/8kuNTDzhDsYwcOZLu3btz0UUXccYZZ4Qsd8kll7Bq1Sqio6OZMGFClv0XXngh9evX57zzzmPo0KEZOgsMGjSIZs2acf3114d7K3wluT1K5vvEIm2Bx1W1i7v+IICqPhNUZqZbZoGInARsB2oAZwDfqGrYnbxjYmI0NjY23/Fu2HGQf4/8kWFXteDixjXzfR5jsvPdhu8AuKzBZT5H4r3Vq1dz7rnn+hrD4cOHqVChAvHx8bRr146RI0eG7NlVnIT6uYjIElWN8eJ6Xv65XgvYHLQeB7TJroyqpojIAaCau6++iCwFDgLDVPVHD2PlSFL6E4xVkRlvPPnDk0DxSDCFwaBBg1i1ahWJiYkMGDCg2CcXP3iZYEK1dmV+XMquzDagrqruEZGWwNci0kRVD2Y4WGQQMAigbt26xxXskUSnXtOqyIxXxl4x1u8QipVPPvnE7xCKPS8b+eOAOkHrtYGt2ZVxq8gqA3tVNUlV9wCo6hJgPXB25guo6khVjVHVmBo1ahxXsIfdBGO9yIxX6lSuQ53KdXIvaEwR4WWCWQw0FJH6IlIK6AdMzlRmMjDAXe4DzFFVFZEabicBRKQB0BDY4GGsQVVk9gRjvDFj3QxmrJvhdxjGRIxnn6Zum8oQYCYQBYxS1ZUiMhyIVdXJwAfAWBFZB+zFSUIA7YDhIpICpAK3q+per2KF4Coye4Ix3nj2p2cB6HpWV58jMSYyPP1zXVWnAdMybXs0aDkRuDrEcROBiV7GltnhxGRKn1SCklE2uIHxxvg+4/0OwZiIsk9T15HEFHt6MZ46rcJpnFbhNL/DKDZee+01mjZtSpMmTXj11VcD2/fu3UunTp1o2LAhnTp1Yt++fQBMnDiRJk2acPHFF7Nnzx4A1q9fT79+/UKe3ysFMQx/hQoVCiia42MJxnUkyeaCMd6asmYKU9ZM8TuMYmHFihW89957LFq0iN9++41vvvmGtWvXAs6gkR07dmTt2rV07NiRZ591qi5feuklFi5cyI033hjogTZs2DBGjBgR1jVTU1O9+WZOYJZgXIcTU6yB33jqpQUv8dKCl/wOo1hYvXo1559/PuXKleOkk06iffv2fPXVV4Az7MuAAU7fogEDBgRGTi5RogRJSUnEx8dTsmRJfvzxR2rWrEnDhg2zvU6FChV49NFHadOmDQsWLGDJkiW0b9+eli1b0qVLF7Zt2wbAe++9R6tWrWjevDlXXXUV8fHxAOzYsYMrrriC5s2b07x5c+bPdyakSx+Gv0mTJnTu3Dkwxtj69evp2rUrLVu25OKLL+aPP/4AYOPGjbRt25ZWrVrxyCOPeHBH80lVi8RXy5Yt9XgMee9HfWjcL8d1DmNysuvILt11ZJffYUTEqlWrMqy3H91eRy8draqqR1OOavvR7XXsb2NVVfXI0SPafnR7Hf/7eFVV3Z+wX9uPbq8TV01UVee+tR/dXif/MVlVVbcd2hbW9Rs2bKi7d+/WI0eO6Pnnn69DhgxRVdXKlStnKFulShVVVf3222+1RYsW2qNHD92/f7927txZ9+7dm+N1AJ0wYYLzfR09qm3bttWdO3eqqur48eP1pptuUlXV3bt3B455+OGH9fXXX1dV1WuuuUZfeeUVVVVNSUnR/fv368aNGzUqKkqXLl2qqqpXX321jh3r3KtLL71U//zzT1VVXbhwoV5yySWqqtqzZ08dM2aMqqq+8cYbWr58+WzvS4jvIVY9+ly2P9ldh5OSOf3knMcaMuZ4VC9X3e8Qio1zzz2X+++/n06dOlGhQgWaN2/OSSfl/HHXqVMnOnXqBMCYMWPo1q0ba9as4cUXX6Rq1aq89tprlCtXLsMxUVFRXHXVVQCsWbOGFStWBM6RmppKzZrOsFMrVqxg2LBh7N+/n8OHD9OlSxcA5syZw0cffRQ4V+XKldm3b1/IYfgPHz7M/PnzufrqY/2ikpKSAPj555+ZONHpF9W/f3/uv//+/N+8AmQJxnXEqsiMx75c/SUAV557pc+RRF7wFAUlo0pmWC9XslyG9cplKmdYr16ueob1cDtK3HLLLdxyyy0APPTQQ9SuXRuAU089lW3btlGzZk22bdvGKaeckuG4+Ph4xowZw8yZM+ncuTOTJk3ik08+Ydy4cdx6660ZypYpU4aoqCjAqQ1q0qQJCxYsyBLLwIED+frrr2nevDkffvghc3OZfiDUMPxpaWlUqVIlMOFYZvmdKsBL1gaD84txJDHZ3uI3nnr9l9d5/ZfX/Q6j2Ni5cycAf//9N19++SXXXnst4EzuNWbMGMB5Uundu3eG455//nnuuusuSpYsSUJCAiJCiRIlAu0m2TnnnHPYtWtXIMEkJyezcuVKAA4dOkTNmjVJTk5m3LhxgWM6duzI22+/DThPPAcPHsx6YlelSpWoX78+n3/+OeB8bv3222+AMwJz8JQDhYUlGCApJY2UNLUnGOOpSf0mManfJL/DKDauuuoqGjduTM+ePXnzzTepWrUqAA888ACzZs2iYcOGzJo1iwceODZV1datW4mNjQ0knXvvvZfzzz+fMWPGBGa9zE6pUqX44osvuP/++2nevDnR0dGBRvsRI0bQpk0bOnXqRKNGjQLHvPbaa3z//fecd955tGzZMpCQsjNu3Dg++OADmjdvTpMmTZg0aVLgPG+++SatWrXiwIEDeb9ZHvFsuP5IO57h+vccSuS6V2fzn25N6dEy9BwOxpjwFYbh+k1WkR6u355gODZMjFWRGS9NWDGBCSuyTjBlTFFldULAYRvo0kTA27FOXXvfpn19jsSYyLBPVGygSxMZ066flnshY4oQSzA4XZQBype222G8U65kudwLFSGqWii7zhZXfrS3WxsMkJDsJJiypSzBGO98vPxjPl7+sd9hRESZMmXYs2ePLx9qJitVZc+ePZQpUyai17VPVCDxqDNIXZlSUT5HYoqy9399H4Abmt3gcyTeq127NnFxcezatcvvUIyrTJkygZdNI8XTBCMiXYHXcCYce19Vn820vzTwEdAS2AP0VdVN7r4HgVtwJhy7U1VnehVnwlHnCaZMSUswxjuz+s/yO4SIKVmyJPXr1/c7DOMzz6rI3CmP3wQuBxoD14pI40zFbgH2qepZwCvAc+6xjXFmt2wCdAXeSp9C2QuJyalElRCbbMx4qmRUSUpGWUcSU3x4+YnaGlinqhtU9SgwHuidqUxvYIy7/AXQUZxWwd7AeFVNUtWNwDr3fNk6nqrexKOplCkZZQ2SxlMfLvuQD5d96HcYxkSMl1VktYDNQetxQJvsyqhqiogcAKq52xdmOrZWThfbsOMgb0xfQZXypVFVWp1Vg0a1qoYVaGJyirW/GM+lJ5eB0QN9jcOYSPEywYR6HMj8nJFdmXCORUQGAYPc1aT/dDtvRZ4izOTT/x7P0WGrDuyOyJWOj8VZsAJxyk2F+kn5hLufhdiJECPAOV6d2MsEEwfUCVqvDWzNpkyciJwEVAb2hnksqjoSGAkgIrFejadTkCzOgmVxFiyLs+CcCDGCE6dX5/ayDWYx0FBE6otIKZxG+8mZykwGBrjLfYA57gxrk4F+IlJaROoDDYFFHsZqjDGmgHn2BOO2qQwBZuJ0Ux6lqitFZDjOFJ2TgQ+AsSKyDufJpZ977EoR+QxYBaQAd6hqqlexGmOMKXievgejqtOAaZm2PRq0nAhcnfk4d99TwFN5uNzI/MToA4uzYFmcBcviLDgnQozgYZxFZj4YY4wxhYu9WWiMMcYTRSLBiEhXEVkjIutE5IHcj/Akhk0i8ruILEvvlSEiJ4vILBFZ6/5b1d0uIvK6G+9yEWkRdJ4Bbvm1IjIgu+vlIa5RIrJTRFYEbSuwuESkpft9r3OPzXMf3GxifFxEtrj3c5mIdAva96B7vTUi0iVoe8jfA7ejyS9u7BPcTid5JiJ1ROR7EVktIitF5C53e2G7n9nFWajuqYiUEZFFIvKbG+cTOZ1bnE4/E9xYfhGRevmNvwBi/FBENgbdy2h3uy8/86BzRYnIUhH5xl33916q6gn9hdOBYD3QACgF/AY09iGOTUD1TNueBx5wlx8AnnOXuwHTcd73OR/4xd1+MrDB/bequ1z1OONqB7QAVngRF07vvrbuMdOBywsoxseBoSHKNnZ/xqWB+u7PPiqn3wPgM6Cfu/wO8O983suaQAt3uSLwpxtPYbuf2cVZqO6p+z1WcJdLAr+49ynkuYHBwDvucj9gQn7jL4AYPwT6hCjvy8886Pr3AJ8A3+T0c4rUvSwKTzDhDEnjl+ChcMYA/wza/pE6FgJVRKQm0AWYpap7VXUfMAtnLLZ8U9UfcHroFXhc7r5KqrpAnd/Oj4LOdbwxZie7YYRC/h64fw1eijMUUebvN69xblPVX93lQ8BqnBEmCtv9zC7O7PhyT937cthdLel+aQ7nzuvQUsf92ZBDjNnx5WcOICK1ge7A++56Tj+niNzLopBgQg1Jk+OwMh5R4FsRWSLOCAMAp6rqNnD+0wOnuNuzizlS30tBxVXLXfYq3iFuNcMocaud8hFjNWC/qqYUZIxulcI/cP6iLbT3M1OcUMjuqVulswzYifOhuz6Hc2cYWgoIHlrKs/9PmWNU1fR7+ZR7L18RZ2T4DDGGGUtB/sxfBf4PSHPXc/o5ReReFoUEE9awMhFwoaq2wBk9+g4RaZdD2eMaIsdDeY3Ly3jfBs4EooFtwEvudt9jFJEKwETgblU9mFPRPMZUoLGGiLPQ3VNVTVXVaJzROloD5+Zwbl/izByjiDQFHgQaAa1wqr3u9zNGEekB7FTVJcGbczh3ROIsCgkmrGFlvKaqW91/dwJf4fxn2eE+AuP+u9Mtnl3MkfpeCiquOHe5wONV1R3uf+w04D2Ojaad1xh341RTnJRpe76ISEmcD+1xqvqlu7nQ3c9QcRbWe+rGth+Yi9Nukd25A/FIeENLFej/p6AYu7rVkKqqScBo8n8vC+pnfiHQS0Q24VRfXYrzROPvvcytkaawf+G8LLoBp0EqvfGpSYRjKA9UDFqej9N28gIZG3+fd5e7k7EhcJEeawjciNMIWNVdPrkA4qtHxgb0AosLZ0ig8znWQNmtgGKsGbT8X5x6YXDmCApuhNyA0wCZ7e8B8DkZGzoH5zNGwakjfzXT9kJ1P3OIs1DdU6AGUMVdLgv8CPTI7tzAHWRsmP4sv/EXQIw1g+71q8Czfv8fCoq5A8ca+X29lxH7EPbyC6fnxp849bcP+3D9Bu4N/w1YmR4DTp3mbGCt+2/6L5TgTMa2HvgdiAk61804DWvrgJsKILZPcapDknH+CrmlIOMCYoAV7jFv4L68WwAxjnVjWI4zNl3wh+PD7vXWENTjJrvfA/fns8iN/XOgdD7v5UU41QLLgWXuV7dCeD+zi7NQ3VOgGbDUjWcF8GhO5wbKuOvr3P0N8ht/AcQ4x72XK4CPOdbTzJefeaaYO3Aswfh6L+1NfmOMMZ4oCm0wxhhjCiFLMMYYYzxhCcYYY4wnLMEYY4zxhCUYY4wxnrAEY4olEZkrIp7Ply4id4ozqvG4TNtjROR1d7mDiFxQgNesJyLXhbqWMZHk6YyWxhRFInKSHhvfKTeDcd4l2Bi8UVVjgVh3tQNwGOcF3YKIoR5wHc6oupmvZUzE2BOMKbTcv8RXi8h77lwc34pIWXdf4AlERKq7Q2QgIgNF5GsRmeLO1zFERO5x58hYKCInB13iBhGZLyIrRKS1e3x5dyDIxe4xvYPO+7mITAG+DRHrPe55VojI3e62d3BedJssIv/NVL6DiHzjDkZ5O/BfceYVuVhEaojIRDeGxSJyoXvM4yIyUkS+BT5y78+PIvKr+5X+FPQscLF7vv+mX8s9x8nu/Vnu3o9mQece5d7XDSJyZ9D9mCrOfCgrRKTv8f1UTbFyPG+M2pd9efmF85d4ChDtrn8G3OAuz8V9SxqoDmxylwfivJ1cEWeYjwPA7e6+V3AGfkw//j13uR3uMDXA00HXqILz5nJ597xxhBi6B2iJ89Z2eaACzmgO/3D3bSLTPEHu9g4ce9v6cYLmacF58sxjPN0AAAJaSURBVLjIXa4LrA4qtwQo666XA8q4yw2B2MznDnGt/wGPucuXAsuCzj0fZ4iQ6sAenKHpr0q/T265yn7/XtjXifNlVWSmsNuoqsvc5SU4SSc336szD8ohETkATHG3/44z9Ee6T8GZj0ZEKolIFaAzzqCBQ90yZXA+5MGdzyPE9S4CvlLVIwAi8iVwMc4QI/lxGdBYjk1sWElEKrrLk1U1wV0uCbwhzmyKqcDZYZz7IpykgarOEZFqIlLZ3TdVncEbk0RkJ3Aqzj17UUSew0lSP+bzezLFkCUYU9glBS2n4gw4CM6TTXoVb5kcjkkLWk8j4+985nGS0oclv0pV1wTvEJE2wJFsYsz3FLfZKAG0DUok6TGQKYb/AjuA5u4xiWGcO6dh1zPf65NU9U8RaYkzDtUzIvKtqg4P67swxZ61wZgT1SacqimAPvk8R18AEbkIOKCqB4CZwH/E/TQXkX+EcZ4fgH+KSDkRKQ9cgTPqbrgO4VTppfsWGJK+4j6hhFIZ2KbO8Pv9cUa9DXW+zLFe7563A7Bbc5jTRkROB+JV9WPgRZyprY0JiyUYc6J6Efi3iMzHaTPIj33u8e/gjOAMMAKn6mm5iKxw13OkzvTEH+KMSvsL8L6q5qV6bApwRXojP3AnEOM2xK/C6QQQylvAABFZiFM9lv50sxxIcRvm/5vpmMfTz43TGWBALrGdBywSZ0bHh4En8/B9mWLORlM2xhjjCXuCMcYY4wlLMMYYYzxhCcYYY4wnLMEYY4zxhCUYY4wxnrAEY4wxxhOWYIwxxnjCEowxxhhP/D9Ei8RPRVCFngAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.axhline(mi,label='ground truth',linestyle='--',color='red')\n",
    "for i in range(rep):\n",
    "    plt.plot(smooth_mi_list[i,:],color='steelblue')\n",
    "    for t in range(smooth_mi_list[i].shape[0]):\n",
    "        if (smooth_mi_list[0,t]>.9*mi):\n",
    "            plt.axvline(t,label='90% reached',linestyle=':',color='green')\n",
    "            break\n",
    "#plt.title(\"Plot of MI estimates against number of iteractions\")\n",
    "plt.xlim((0,smooth_mi_list[0].shape[0]))\n",
    "plt.ylim((0,mi*1.1))\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"MI estimate\")\n",
    "plt.legend()\n",
    "plt.savefig(fig_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
